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ABSTRACT 


At',  automated  general  purpose  system  for  analysis  is 
presented.  This  system,  identified  by  the  acronym  "MAGIC  TI" 
for  Matrix  Analysis  via  Generative  and  Interpretive  Computations, 
is  an  extension  of  structural  analysis  capability  available  in 
the  initial  MAGIC  System.  MAGIC  provides  a  powerful  framework 
for  implementation  of  the  finite  element  analysis  technology  and 
provides  diversified  capability  for  displacement,  stress,  vibration 
and  stability  analyses. 

The  matrix  displacement  method  of  analysis  based  upon 
finite  element  idealization  is  employed  throughout.  Ten  versatile 
finite  elements  are  incorporated  in  the  finite  element  library. 

These  are  frame,  shear  panel,  triangular  cross-section  ring, 
trapezoidal  cross-section  ring  (and  core),  toroidal  thin  shell 
ring  (and  shell  cap),  quadrilateral  thin  shell  and  triangular  thin 
shell  elements.  Additional  elements  include  a  frame  element, 
quadrilateral  plate  and  triangular  plate  elements  which  can  be 
used  for  both  stress  and  stability  analysis.  The  finite  elements 
listed  include  matrices  for  stiffness,  mass,  incremental  stiffness, 
prestrain  load,  thermal  load,  distributed  mechanical  load  and  stress. 

The  MAGIC  II  System  for  structural  analysis  is  presented 
as  an  integral  part  of  the  overall  design  cycle.  Considerations  in 
this  regard  include,  among  other  things,  preprinted  input  data  forms, 
automated  data  generation,  data  confirmation  features,  restart _ 
options,  automated  output  data  reduction  and  readable  output  displays. 

Documentation  of  the  MAGIC  I”  System  is  presented  in 
three  parts;  namely.  Volume  I:  Engin  er's  Manual,  Volume  II:  User's 
Manual  and  Volume  III:  Programmer's  Manual.  The  subject  document. 
Volume  I  (Engineer's  Manua]  -  Addendum)  is  an  extension  of  the 
primary  technical  document.  Included  are  the  theoretical  develop¬ 
ments  for  the  additional  finite  elements  included  in  the  MAGIC  II 
System  as  well  as  a  discussion  of  newly  added  computational  procedures. 
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SECTION  I 


INTRODUCTION 

'T’he  MAGIC  II  Systems  for  Structural  Analysis  is  a  logical 
ext.jrition  of  the  original  MAGIC  System  reported  in  References  1, 

?  and  3.  All  capabilities  available  from  the  original  MAGIC  System 
nave  been  letained.  Extension  of  the  program  capability  is  primarily 
in  the  following  areas. 

(a)  The  implementation  of  four  additional  finite  element 
representations  and  their  associated  element  matrices. 

(b)  The  improvement  of  output  displays  to  facilitate 
ease  of  interpretation  by  the  User. 

(c)  The  provision  of  an  "Agendum  Library"  to  accommodate 
the  following  classes  of  analyses. 

(1)  Statics 

(2)  Statics  with  Condensation 

(3)  Statics  with  Prescribed  Displacements 

(4)  Stability 

(5)  Dynamics  (Modes  and  Frequencies) 

(6)  Dynamics  (with  Condensation) 

(d)  The  addition  of  an  out-of-core  eigenvalue  routine 

for  the  nonsymraetric  matrices  based  on  the  power  method 
"on  the  order  of”  3000  x  3000. 

(e)  The  addition  of  Improved  and  expanded  error  diagnostics. 

(f)  The  addition  of  a  prescribed  displacement  option  to 
accommodate  more  than  one  load  condition  per  execution. 

(g)  The  addition  of  the  capability  to  accept  either 
rectangular,  cylindrical  or  spherical  coordinates  as 
input  data. 

(h)  The  addition  of  miscellaneous  arithmetic  modules  to 
the  System  to  support  the  computational  procedures. 
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(i)  The  addition  of  a  new  assembly  module  to  increase 
the  permissible  assembled  system  matrix  size. 

Documentation  of  the  ['1AGIC  II  System  for  structural  analysis 
is  presented  in  three  volumes.  The  subject  volume  (Volume  I)  is  an 
addendum  to  the  primary  technical  report  documented  in  Reference  1. 
Separate  supplementary  volumes  are  provided  to  facilitate 
utilization  of  the  MAGIC  II  System.  Volume  II,  the  User's  Manual^^^, 
Includes  detailed  specifications  for  the  preparation  of  input  data, 
along  vjith  illustrative  examples.  Volume  III,  the  Programmer's 
Manual  presents  Information  on  the  organization  of  the  MAGIC  II 

System  as  well  as  its  operational  characteristics. 

It  is  to  be  noted  that  this  addendum  is  to  be  used  in 
conjunction  with  the  original  technical  report  (Reference  1)  in 
order  to  utilize  the  MAGIC  II  System  effectively. 

Section  II  presents  a  discussion  of  additional  analysis 
and  programming  technology  which  has  been  incorporated  into  the 
MAGIC  II  System.  New  computational  procedures  and  expanded 
size  characteristics  are  emphasized. 

A  general  theoretical  description  of  the  additional 
finite  element  representations  (and  element  matrices)  included  in 
the  MAGIC  II  System  is  given  in  Section  III.  These  elements  are 
as  follows ; 

(a)  Trapezoidal  Cross-Section  Ring  (Core) 

(b)  Quadrilateral  Plate 

(c)  Triangular  Plate 

(d)  Incremental  Frame 

Section  IV  presents  a  discussion  of  the  computational 
procedures  available  in  the  MAGIC  II  System.  Included  are  procedures 
for  Statics,  Statics  with  Condensation,  Statics  with  Prescribed 
Displacements,  Stability,  Dynamics  (Modes  and  Frequencies)  and 
Dynamics  with  Condensation.  Additional  procedures  are  outlined  for 
Static  and  Dynamic  Substructuring. 
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"'ne,  body  oC  this  technical  report  is  concluded  with  a 
general  retrospective  discussion  in  Section  V.  Limitations  of 
t,hi-  MAiliC  1  I  System  are  discussed  and  guidelines  for  utilization 
are  presented. 
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SECTION  II 


ADDITIONAL  ANALYSIS  AND  PROGRAMMING  TECHNOLOGY 

A.  ANALYSIS  TECHNOLOGY 

The  MAGIC  IF  System  incorporates  the  ten  finite  elements 
shown  in  Figures  II.l  and  I I. 2  .  The  six  finite  elements  shown 

in  Figure  II.l  were  available  in  the  original  MAGIC  System  and 
are  discussed  in  detail  in  Reference  1.  The  four  additional  finite 
elements  shown  in  Figure  II. 2  are  described  in  detail  in  Section  IV. 

The  set  of  matrices  embodied  in  each  element  representation 
determines  the  type  of  analyses  which  can  be  performed.  In  the 
MAGIC  II  System,  a  complete  element  representation  is  taken  to 
include  matrices  for  stiffness,  incremental  stiffness,  pressure 
load,  prestrain  load,  thermal  load,  stress,  and  mass.  Moreover, 
provision  has  been  made  for  additional  element  matrices  such  as 
consistent  damping  matrices. 

The  types  of  analyses  available  with  the  MAGIC  II  System 
are  as  follows; 

(a)  Statics 

(b)  Statics  (With  Condensation) 

(c)  Statics  With  Prescribed  Displacements 

(d)  Stability 

(e)  Dynamics  (Modes  and  Frequencies) 

(f)  Dynamics  (With  Condensation) 

In  addition  many  user  variations  of  the  above  computational 
procedures  are  available  with  the  System.  This  is  possible  due 
to  the  powerful  matrix  abstraction  capability  available  from  the 
MAGIC  II  System.  A  complete  description  of  the  computational 
procedures  listed  above  along  with  example  problems  which  demonstrate 
their  use  is  provided  in  Volume  II  of  this  report  (User's  Manual). 
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Figure  HJMAGIC  System  Finite  Elements 


Triangular  Plate  Element 


B.  PROGRAMMING  TECHNOLOGY 

In  this  section  additional  programming  technology  available 
with  t’ne  MAGIC  II  System  is  discussed.  Volume  III  of  this 
report  (Programmer's  Manual)  is  suggested  for  complete  documenta¬ 
tion  on  program  technology. 

r.  General  Description 

The  general  arrangement  of  the  MAGIC  II  digital  computer 
program  system  is  shown  in  Figure  n.  3  .  The  supervisory  program 
consists  of  the  FORMAT  control  and  two  monitors;  the  Preprocessor 
Monitor,  and  the  Execution  Monitor.  The  main  program  controls 
the  normal  two  phase  operation  by  delegating  control,  in  turn, 
to  the  two  monitors. 

The  preprocessor  Monitor  directs  the  processing  of  card 
input  data  describing  the  machine  configuration,  the  problem 
specification,  the  abstraction  instruction  sequence  and  the  matrix 
data. 

A  standard,  modified  standard,  or  totally  new  machine 
configuration  may  be  defined  for  each  MAGIC  II  case. 

General  output  format  and  labeling  information,  and  identify¬ 
ing  names  of  the  master  input  and  output  data  sets  (tapes)  con¬ 
stitute  the  problem  specification  data. 

The  matrix  and  pseudo-matrix  operations  are  input  in  the 
required  sequence  of  execution  in  the  abstraction  instruction 
sequence.  Abstraction  instructions  are  submitted  in  free  form  on 
standard  Fortran  coding  sheets  for  punched  card  reproduction. 

Card  input  matrix  data  are  specified  on  a  standard  form. 
Matrices  may  be  of  order  3000x3000,  and  may  contain  up  to  6OOO 
randomly  ordered,  single  precision  real  elements  (with  32K  core 
storage  unit). 
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Figure  II. 3  MAGIC  II  -  Digital  Computer  Program 


For  the  general  case,  preprocessing  involves  straightforward 
sequential  processing  of  data  by  each  of  the  modules  under  the 
Preprocessor  Monitor.  Special  preprocessing  can  be  specified  by 
proper  use  of  the  control  cards  described  in  the  User's  Manual. 

The  final  preprocessor  operation  is  to  pre-plan  the  data 
storage  allocation  through  the  problem  and  to  record  this  program 
of  the  "complete  problem  solution  logic"  for  use  by  the  Execution 
Monitor. 

The  standard  matrix  operational  modules  provide  for  matrix 
addition,  subtraction,  multiplication,  and  transpose  multiplication, 
with  optional  concurrent  scaling,  and  for  matrix  scalar  multiplication, 
transposition,  adjoining,  de joining,  and  inversion.  Modules 
for  the  solution  of  simultaneous  equations  by  elimination  and  iterative 
techniques  complete  the  basic  standard  matrix  operation  capability 
of  the  system. 

The  pseudo-matrix  operational  modules  provide  for  the  element 
by  element  multiplication  of  two  matrices  of  Identical  order,  the 
elements  of  a  matrix  to  be  raised  to  a  scalar  power,  the  extraction 
of  the  algebraic  maximum  and  minimum  elements  of  the  rows  or  columns 
of  a  matrix  (i.e.,  the  envelope  of  a  matrix),  the  diagonal izat ion  of 
a  row  or  column  matrix,  the  generating  of  null  and  identity  matrices, 
and  the  renaming  of  a  matrix.  Included  in  the  classification  of 
pseudo-matrix  operational  modules  is  the  "Structure  Cutter"  subroutine 
which  generates  a  well  conditioned  solution  of  "n"  linear 
simultaneous  equations  in  "m"  unknowns  by  Jordaniein  elimination 
(where  n  i  m) . 

Matrices  produced  as  the  results  of  standard  and  pseudo-matrix 
operations  may  be  as  large  as  3000x3000  with  no  restriction  on 
population  density.  Storage  of  matrix  data  is  by  column  sort,  and 
when  individual  column  population  density  is  less  than  50  percent, 
storage  is  in  compressed  format.  In  compressed  format,  each  non¬ 
zero  element  and  its  corresponding  row  location  are  sequentially 
stored,  and  zero  elements  are  omitted.  Where  feasible,  the  sub¬ 
routines  operate  directly  on  the  compressed  data. 
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MAGIC  IT  includes  two  subroutines  for  the  calculation  of 
eigenvalues.  The  first  subroutine  calculates  the  specified  number 
of  eigenvalues,  beginning  with  the  largest,  and  the  corresponding 
eigenvalues  of  a  matrix,  whose  maximum  order  is  limited  by  the 
workiiig  core  storage  availaole  to  the  subroutine.  Typically,  with 
a  3?K  storage  unit,  the  matrix  may  be  as  large  as  l60  x  l60. 

This  subroutine  is  written  for  a  real  symmetric  matrices  only.  The 
second  subroutine  also  calculates  the  specified  number  of  eigenvalues 
and  eigenvectors  beginning  with  the  largest  eigenvector.  However, 
the  real  eigenm.atrix  can  be  symmetric  or  nonsymmetrlc  and  the  only 
limit  on  its  order  is  the  amount  of  working  storage  available  to  the 
MAGIC  System. 


Up  to  nine  special  operational  subroutines  can  be  coded  by  the 
user  and  added  to  the  system.  The  fourth  user  coded  module  is  the 
structural  generative  system  of  MAGIC  and  is  described  in  detail 
in  the  User’s  Manual. 

The  sequence  of  operation  is  controlled  by  simple  abstraction 
instructions  prepared  by  the  User,  keypunched,  and  read  directly 
by  the  machine.  Comments  may  be  included  in  the  abstraction  in¬ 
struction  sequence  for  explanation  of  the  results. 

Limited  logic  is  available  in  the  form  of  a  conditional  trans¬ 
fer.  A  matrix  may  be  tested  for  nullity  and,  if  true,  control  will 
be  transferred  forward  to  a  specified  abstraction  instruction  in 
the  sequence.  Conditional  transfer  is  limited  to  a  "skip  ahead" 
in  the  abstraction  instruction  sequence. 

Matrices  can  be  printed  in  a  standard  form,  with  small  number 
suppression  and  row-column  labeling.  The  matrix  elements  are 
printed  as  floating  point  numbers  with  optional  exponent. 

The  normal  printed  output  for  a  MAGIC  II  case  includes  a 
listing  produced  by  the  preprocessor.  The  listing  unconditionally 
includes  all  control  and  specification  data  together  with  the  complet 
abstraction  instruction  sequence.  The  listing  will  also  include 
matrix  input  data,  special  input  data,  and  the  machine  generated 
"complete  problem  solution  logic"  if  the  approriate  options  are 
chosen  in  the  control  data. 


II. 


structural  Abstraction  Instructions 


In  designing  the  MAGIC  II  System  for  Structural  Analysis, 
provision  was  made  for  accommodating  new  abstraction  instructions 
peculiar  to  the  .USER04.  module.  In  keeping  with  the  philosophy  of 
generating  a  highly  flexible  USER  oriented  system,  specialized 
instructions  were  designed  for  items  such  as  element  stress  and 
force  determination,  element  assembly  and  print  controls.  These 
additional  USER  options  provide  output  capabilities  of  the  MAGIC  II 
System,  consistent  with  input  requirements. 

The  following  abstraction  instructions,  .STRESS.,  .FORCE., 
.ASSEM.,  .EPRINT.,  and  .GPRINT.  are  to  be  used  in  conjunction  with 
the  .USER04.  abstraction  Instruction. 

To  compute  the  net  element  stress  matrix  and  generate 
optional  engineering  print  of  apparent  element  stresses,  element 
applied  stresses  and  net  element  stresses  use  the  .STRESS, 
abstraction  instruction. 

To  compute  the  net  element  force  matrix  and  generate 
optional  engineering  print  of  apparent  element  forces,  element 
applied  forces  and  net  element  forces  use  the  .FORCE,  abstraction 
instruction. 

To  assemble  the  element  stiffness  matrices,  element  mass 
matrices,  element  incremental  matrices  and  element  thermal  load 
matrices  as  output  by  the  .USER04.  Instruction  use  the  .ASSEM. 
abstraction  instruction. 

To  generate  engineering  printout  of  the  net  element  stresses 
or  net  element  forces  use  the  .EPRINT.  abstraction  instruction. 

To  generate  engineering  printout  of  reactions,  displacements, 
eigenvalues  and  eigenvectors,  and  user  matrices  use  the  .GPRINT. 
abstraction  instruction. 

A  complete  discussion  of  the  above  listed  instructions  along 
with  a  detailed  explanation  of  their  proper  usage  is  presented 
in  Section  II.B.f  of  the  User’s  Manual. 
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III.  Size  Characteristics 

The  size  characteristics  of  the  MAGIC  II  System  are 
twofold:  first,  there  are  the  size  characteristics  of  the  program 
itself  and  second,  those  associated  with  the  problem  solving 
capability.  Considering  the  former,  the  MAGIC  II  System  contains 
356  subroutines  (approximately  38OOO  FORTRAN  IV  source  cards) 
logically  designed  into  43  overlay  links  on  an  IBM  36O/65 
using  45600  words  of  storage.  The  overlay  design  reflects  the 
optimum  use  of  available  storage,  yet  maintains  respectable 
execution  efficiency. 

The  MAGIC  II  System  offers  large  scale  capability  with  no 
penalties  to  small  applications  due  to  the  fact  that  out  of  core 
operations  are  not  utilized  unless  the  magnitude  of  the  application 
requires  them. 

The  scale  of  the  analysis  capability  provided  via  the  MAGIC 
II  System  can  be  characterized  as  "on  the  order  of"  3000  displacement 
degrees -of -freedom  using  45600  words  of  storage  on  an  IBM  36O/65 
digital  computer.  Other  relevant  maximum  size  characteristics 
are  3000  discrete  elements  and  1000  gridpoints.  Matrices  which 
are  card  input  may  be  of  order  3000  x  3000  and  contain  up  to 
6000  single-precision  real  non-zero  elements. 

The  MAGIC  II  System  needs  a  minimum  of  eight  external  storage 
units  to  operate,  distributed  into  the  following  functions:  one  unit 
assigned  as  Instruction  storage  for  the  Execution  Monitor,  one  unit 
assigned  as  a  Master  Input  Unit,  one  unit  assigned  as  a  Master 
Output  Unit,  and  five  units  assigned  as  Input/Output  Utility  Units. 
Every  effort  should  be  made  to  make  the  most  external  storage  units 
possible  available,  since  any  increase  in  the  available  storage  units 
Increases  execution  efficiency. 
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SECTION  III 


ADDITIONAL  FINITE  ELEMENT  REPRESENTATIONS 

A.  INTRf^DUCTION 

The  MAGIC  II  System  incorporates  ten  finite  element 
representations.  Six  of  these  elements,  namely,  frame,  shear 
panel,  triangular  cross-section  ring,  toroidal  thin  shell  ring, 
quadrilateral  thin  shell  and  triangular  thin  shell  were 
available  in  the  initial  MAGIC  System  and  are  described  in 
detail  in  Reference  1. 

Four  additional  elements,  namely,  trapezoidal  cross-section 
ring  (core),  quadrilateral  plate,  triangular  plate  and  incremental 
frame  have  been  incorporated  into  MAGIC  II.  A  complete  element 
representation  is  taken  to  include  matrices  for  stiffness, 
stress,  incremental  stiffness,  pressure  load,  prestrain  load, 
thermal  load  and  mass. 

In  the  following  sections,  each  of  the  element  representations 
along  with  associated  element  matrices  are  discussed  in  detail. 


13 


B.  TRAPEZOIDAL  CROSS-SECTION  RING  (CORE) 

I.  Introduction 

The  formulation  of  the  trapezoidal  cross-section  ring 
discrete  element  described  herein,  if  dervied  from,  and  is 
mathematically  consistent  with,  the  formulation  described 
in  Reference  6. 

The  trapezoidal  cross-section  ring  discrete  element,  shown 
in  Fig. 1  r 1 . 1,  provides  a  powerful  tool  for  the  analysis  of  thick 
walled  and  solid  axisymmetric  structures  of  finite  length.  It 
may  be  used  alone  or  if  the  problem  dictates  a  highly  irregular 
grid  work,  it  may  be  combined  with  the  well  known  triangular  ring 
discrete  element  to  form  the  assembly  of  any  axisymmetric 

structure  taking  into  account: 

1.  arbitrary  variations  in  geometry 

2.  axial  variation  in  orientation  of  material  axes 
of  orthotropy 

3.  radial  and  axial  variations  in  material  properties 

4.  any  axisymmetric  loading  system  which  can  include 
pressure,  and  temperature,  and  degradation  of 
material  properties  due  to  temperature. 

For  the  analysis  of  solid  structures,  a  core  discrete  element 
has  been  developed  to  be  used  in  conjunction  with  the  trapezoidal 
or  triangular  cross-section  discrete  element.  This  core  element 
(Fig .III .2)is  a  specilizatlon  of  the  trapezoidal  cross-section 
ring  element. 

The  discrete  element  technique  was  first  applied  to  the 
analysis  of  axisymmetric  solids  by  Clough  and  Rashid^S^)  and 
later  the  formulation  of  the  traingular  cross-section  ring  was 
extended  by  Wilson to  include  non-axi symmetric  loading.  This  develof 
ment  deals  with  the  axisymmetric  case  but  includes  orthotropic 
material  properties.  The  integration  of  the  strain-energy  over 
the  volume  of  the  ring  is  effected  analytically,  and  finally 
pre-strain,  and  pressure  load  vectors  as  well  as  a  consistent  mass 
matrix  are  included.  Thus  the  following  element  representation  is 
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Special  Conditions  On  Core  Element 


(a)  n  =  Th  s  0 

(b)  Ul  ^  U4  H  0 


FIGURE  III. 2  CORE  ELEMENT  SPECIALIZATION 
OP  TRAPEZOIDAL  CROSS-SECTION 
RING  DISCRETE  ELEMENT 


•-  v**' 


formulated  to  include  algebraic  expressions  for  the  following 
matrices  : 


1. 

Stiffness 

% 

[K] 

2. 

Stress 

> 

[S] 

3. 

Consistent  Mass 

> 

[M] 

4. 

Pressure  Load 

$ 

{Pp} 

>  • 

Thermal  Load 

> 

{Prp} 

6. 

Gravity  Load 

> 

(Fq) 

7. 

Centrifugal  Force 

> 

The  above  matrices  arise  as  coefficient  matrices  in  the 
generalized  form  of  the  Lagrange  Equations  for  the  element.  The 
generalized  form  of  the  Lagrange  equation  appropriate  for  the 
complete  element  representation  listed  above  is  given  by. 


d 

^  dt 


0 


Qj,  *  generalized  displacement 

q^  -  generalized  velocity 

=  total  potential  energy 
4>2  =  kinetic  energy 
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I  1  Assumed  Displacement  Functions 


A  structural  element  Is  mathematically  discretized  into  a 
finite  number  of  degrees  of  freedom  by  the  assumption  of 
displacement  function  mode  shapes.  The  displacement  modes 
employed  for  the  trapezoidal  ring  may  be  written  as: 


u(r,0,z)  = 

oti 

+  aar  + 

asz  + 

a4rz 

(2.1) 

w(r,e,z)  = 

3i 

+  32r  + 

33Z  + 

34i*z 

(2.2) 

It  is  to  be  noted  that  the  assumed  displacement  functions  are 
interelement  continuous  when  the  elements  employed  are 
rectangular.  The  coefficients  a  and  3  which  appear  in  the 
assumed  displacement  functions  will  be  referred  to  as  field 
coordinate  displacement  degrees  of  freedom.  The  transformation 
from  field  coordinates  to  grid  point  displacement  degrees  of 
freedom  (u^)  is  effected  by  writing: 


a  +  a  r.  + 
1  2  i 


a  z. 
3  i 


a  r.  z. 
4  1  i 


(2.3) 


Hence 


(a)  =  [h]{u} 


also 


'^i^^i®i^i^  =  + 


3  r.  + 

2  1 


3  z,  +  3  r.z. 
3  1  4  1  1 


Then 


{3}  =  [h]{w} 

Upon  combination  of  (2.4)  and  (2.6)  we  have 

iy]  =  [H]{q} 


(2.4) 


(2.5) 


(2.6) 


(2.7) 


18 


where 

(2.8) 

(2.9) 

A  special  case  arises  when  the  trapezoidal  ring  Is  to  be 
used  as  a  core  element.  For  this  (see  Figure  m. 2) 
ri  =  r I,  ^0  and  ui  h  U4  s  0.  This  causes  the  quantities 
ai  and  as  In  the  assumed  displacement  mode  to  be  equal  to 
zero,  which  causes  the  [H]  matrix  to  be  modified.  This 
modified  matrix  Is  designated  [H]  for  the  core  element 
specialization. 


(y)^  =  Loti,  oi2,  as,  aif,  3i,  Ba,  Bs,  B4j 

T  I 

{q}  =  |_Ul,  Wl,  U2,  W2,  Us,  Ws,  U4,  Wi,  J 
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Potential  lilnergy 


1  I  1 

The  total  potential  energy  is  derived  in  this  section  as  the 
sur:i  of  the  strain  energy  and  external  work  contributions. 

The  procedure  followed  is  exactly  the  same  as  that  detailed 
in  Pveference  1.  The  desired 

form  of  the  potential  energy  is  as  follows: 

U  =  I  J [LeJ[E]{e}  -  LeJ[E]{e.}]dV 
V 


(3. 


TV  i'lletiieiiL  Glatic  Matrices 

^ . 1  Introduction 

In  order  to  effect  the  discretization  of  the  element, 
the  assumed  displacement  modes  must  be  introduced  into  the 
potential  energy  function.  Substitution  of  the  total  potential 
energy  function  into  the  Lagrange  equation  yields  the  element 
matrices  with  respect  to  grid  point  displacement  degrees  of 
freedom.  An  exception  is  the  element  stress  matrix  which  is 
derived  from  the  strain-displacement  and  stress-strain 
relationships . 

4 . 2  Stiffness  Matrix 

The  energy  contribution  to  the  linear  elastic 
stiffness  is 


j  J  LeJ[E]{e}dV  (4.2.1) 

V 

The  strains  can  be  expressed  in  terms  of  the  generalized 


coordinates 

using  Equations  (2. 

1),  (2.2), 

and 

the 

fact 

that  {e}'^ 

u/r, 

w^,  u^+w 

r  J 

Then 

ie)  = 

Cd]{y} 

(4.2.2) 

where 

0 

1 

0  z 

0 

0 

0 

0 

1/r 

1 

1  2 

0 

0 

0 

0 

[Dl  = 

(4.2.3) 

0 

0 

0  0 

0 

0 

1 

r 

0 

0 

1  r 

0 

1 

0 

Z 

Since 

LeJ  = 

LvJCb]" 

(  4 . 2  .  it ; 
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We  have  upon  substitution  into  (4.2.1) 


=  1/2  J  LyJ  [D]'^[E][D]{Y}dV  (4.2.5) 

V 

For  the  elemental  volume  of  the  ring  element  in  cylindrical 
coordinates  we  have: 


dV  =  2-irrdrdz  (4.2.6) 

Substituting  back  into  the  strain  energy  equation  we  can  now 
write : 


# 


=  1/2  2TrLxJjJr[D]  [E][D]drdz{Y} 
r  z 


(4.2.7) 


All  of  the  integrals  in  Equations  (4.2.7)  are  of  the  form 


pq 


=  //= 


r^z^drdz 


(4.2.8) 


r  z 


It  is  now  desirable  to  see  how  the  integration  is  carried  out 
over  the  trapezoidal  cross-section. 


r  ,u 


a  = 


• 

% 

b 

_  rif-ri 

Z4-Z1 

2  If  *^2  1 

r2  -  (rT)22 

• 

3 

d 

_  r3-r2 

Z3-Z2 

Z3  -Z2 

!f' 


r  I 


.  s 


s 


For  the  trapezoid 


I 

pq 


the  integration  takes  the  'form: 


For  the  case  with  the  side  r  =  c  +  dz  parallel  to  the  axis  of 
symmetry  (the  7.  axis)  we  have: 

Zi,  c 

Ipq  =f(  r^z^drdz  (4 

a+bz 


For  the  case  with  the  side  r  =  a  +  bz  parallel  to  the  axis  of 
symmetry  we  have: 


Z>f 

z  1  a 


c+dz 


■Pz^^drdz 


r^'z 


And  finally  for  the  rectangle,  the  integration  takes  the  form: 

c 

I  =rr  r-Pz^drdz 

pq  J  J 

zi  a 


For  the  case  where  r  =  c  +  dz  t  c  and 

r  =  a  +  bz  -*•  a  . 

a  test  is  made  in  the  computer  and  we  have  the  following 
Let 

’  = 

I 
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.2.9) 


.2.10) 


.2.11) 


.2.12) 


I  I 


il 

Then 

By  the  same  token 
Let 


ri-r4 

d 


ri  =  d 


<  e  (prescribed) 
and  r4  =  d 


d* 


r2-tr3 

2 


If 

Then 


r2-r3 

d* 


e 


(prescribed) 


r2  =  d'  and  r3  =  d* 
Equation  (4.2.7)  can  now  be  rewritten  as: 


=  1/2  LYJ  CK]{y} 


(4.2.13) 


where  [K]  is  shown  on  Page  10  of  Reference  6. 

Introducing  the  transformation  to  gridpoint  displacement  degrees 
of  freedom  we  have: 


(y)  = 

[H]{q) 

(2.7) 

LyJ  = 

LqJH]'' 

(4.2.14) 

Then  (4.2.13)  becomes 

=  l/2LqJ[H]'^[K][H]{q}  (4.2.15) 

Upon  taking  the  first  variation  with  respect  to  the  displacements, 
we  obtain  the  element  stiffness  matrix  referenced  to  grid  point 
displacement  degrees  of  freedom 

( 

[K]  =  [H]'^CK][H] 
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(4.2.16) 


For  the  special  case  of  the  core  element,  the 
element  stiffness  matrix  referenced  to  grid  point  displacement 
degrees  of  freedom  is  obtained  as  follows: 

[K]*‘  =  [H]'^CK][H]  (i|.2.17) 

4 . 3  Pressure  Load  Matrix 

The  pressure  load  matrix  will  be  developed  in  the 
following  manner:  The  pressure  load  due  to  pressure  normal  to 
the  sides  between  node  points  (1)  and  (^),  and  (2)  and  (3)  will 
be  developed  first  and  then  the  axial  pressure  load  (the  pressure 
normal  to  the  sides  connecting  nodes  (1)  and  (2)  and  (3)  and  (4) 
will  be  developed  next.  These  will  then  be  combined  so  that 
radial  and  axial  pressures  may  be  input  for  each  node  point  of 
the  trapezoidal  element  (See  Pig  .111.1  for  node  point  numbering). 


^1.3.1  Radial  Pressure 

Assume  a  linear  normal  pressure  distribution  on  the 
boundary  between  node  points  (1)  and  (4).  This  assumption  leads 
to  the  requirement  of  numbering  the  node  points  in  counterclock¬ 
wise  order. 


Let 


where 


p  =  Pi  +  e  +  fz 


e 

f 


4EjlzE4 

(Zh-Zi  ) 


Zl 


EJiZEi 

Z4-Z  1 


(4. 3. 1.1) 


(4. 3. 1.2) 


The  external  work  done  by  the  pressure  on  the 
displacements  is 


W 


I 


w)d.A 


(4. 3. 1.3) 
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(^.3.1.5) 

(i4.3.1.6) 

(^.3.1.7) 

(4. 3. 1.8) 


u(a+bz,z)  =  ai  +  aaCa+bz)  +  asz  +  a^(a+bz)z  (4. 3. 1.9) 

w(a+bz,z)  =  +  B2(a+bz)  +  &3Z  +  3i,(a+bz)z  (^.3.1.10) 

For  the  side  of  the  trapezoid  connecting  node  points  (2)  and 
(3)  the  same  procedure  is  followed  exactly.  The  total  work 
can  be  written  as: 

W  =  ZW  ;  (J  =1,2,  ...,8)  (il.3.1.11) 

S 

W  =  LyJCQ„]{M}  (i4.3.1.12) 

P 
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where 


Recalling  that 

( —  1 — 
3  -<: 

l_  L_ 

II  II 

Lai,a2,a3,a4,6i,B2,33,34j 

l_m j  iHjj  n^jJ 

(4.3.1.13) 

(y)  =  [H]{q} 

(2.7) 

then 

LyJ  =  LqJCh]'^ 

(4.2.14) 

Then 

W  =  LqJ[H]'^[Q  ]{M} 

(4.3.1.14) 

The  vector  {M}  can  be  written  in  the  following  manner: 


/  N 

m_ 

I 

mil 

)  = 

n. 

I 

n__ 

II 

\ 

0  0  -(:r-^) 


Zl 


'  Zi,-Z  1 

0  l+(— ^ 


Z4-Z 1 


■)  -(r-^i-)  0 


zi,-Zi:  'zi,-zj 


-<i;riT> 


0  ( — ± — ) 

'‘Z4-Z1  ^ 


/  \ 
Pi 


P2 

P3 

P- 

V  / 


>(4.3.1.15) 


or 


m  =  [h  ]{p} 

Equation  (4.3.1.14)  can  now  be  written  as 

W  =  LqJCH]'^[Qp][hp]{p} 

where  the  radial  pressure  load  vector  is  {Pp}  and  has  the 
following  form 


{Fp}  =  [H]^[Qp][hp]{p} 


(4.3.1.16) 


(4.3.1.17) 


(4.3.1.18) 
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i4.3.2  Axial  Pressure 

Assume  a  linear  normal  pressure  distribution  on 
the  boundary  between  node  points  (1)  and  (2) 


p.  =  pi  _  (E2z£l)ri  +  (EllEl)r 
^  '‘ra-ri  *  r2-ri 


('1.3.2.1) 


I  t 

P2-P1 

r2“ri 


(^♦.3.2.2) 


Then 


p  =  (mj  +  nj  r) 


(4. 3. 2. 3) 


The  external  work  done  by  the  pressure  on  the 
displacement  is: 


■  / 


W  =  J  (Pj,u  +  p^v)dA 


(4. 3. 2. 4) 


For  the  case  of  axial  pressure;  p^  =  0 
Therefore 


=  J  (PgW] 


(4. 3. 2. 5) 


r^  2tt 
^or 


-n 

r,  0 


rd0dr 


(4. 3. 2. 6) 


2Tr  /  rdr 


A 


(4. 3. 2. 7) 


Substituting  into  Equation  (4. 3. 2. 5)  the  following  result  is 
obtained : 

Ti 

V  =  2v  J  r(p2w)dr 

r  1 

For  the  side  of  the  trapezoid  connecting  node  points  (3)  and 
(A)  the  same  procedure  is  followed  as  for  node  points  (1)  and 
(2).  The  total  work  can  then  be  written  as: 


(i». 


and 


[Q  ] 

Lx^pJ 


W  =  E  W,; 

J  J 

(j  =  1,2, 3, 4) 

W  =  e  [Qp]{M'} 

LBJ  =  L6i,32,B 

3,B4j 

LMJ  =  Uj, 

’  •  1 

l/2(r2-ri)  l/2(r4-r3) 

l/3(r2-r? ) 

l/3(r4-rl) 

l/3(r2-ri)  l/3(r4-r3) 

l/4(r2-ri ) 

l/4(r4-r5) 

2ti 

Zi/2(r2-ri)  z^/2(r^-r3) 

Zi,/3(r2-ri ) 

Z4/3(r2-r3 ) 

Zi/3(r2-r!)  Z4/3(rH-rL 

Zi/4(r2-r![ ) 

Zi,/4(rll-r5) 

known  that 

{&}  = 

Ch]{w} 

LBJ  = 

LwJCh]'^ 

(4. 

(4. 

(4. 


(4, 


(2. 


Substituting  into  Equation  (4.3.2.10)  we  obtain 


W  =  LwJ[h]'^CQp]{K'} 


(4. 


3.2.8) 

3.2.9) 

3.2.10) 

3.2.11) 

3.2.12) 

6) 

3.2.13) 
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f 

The  vector  {M  }  can  be  written  in  the  following  manner: 


f 


t 


s 


t 


r 


ra- 


X 

r 


1 


0 


1 

ra-ri 

0 


r,-ri 


0 


(  \ 

0 

Pi 

r4-r3 

pI 

< 

0 

pI 

1 

Pi* 

1  / 

(4.3.2.14) 


or 

{m'}  =  [hp]{p’} 

Equation  (4.3.2.13)  can  now  be  written  as 


W  =  LwJ  Ch]'^[Qp]Chp]{p’ } 


(4.3.2.15) 


(4.3.2.16) 


{Fp}  =  Ch]^[Qp][hp]{p’} 


(4.3.2.17) 


4.3.3  Combining  Radial  and  Axial  Pressure  Loads 

In  this  section  we  will  combine  the  radial  and  axial 
pressure  load  vectors  so  that  it  is  possible  to  input  one  radial 
and  one  axial  pressure  value  for  each  node  point  of  an  element. 

From  Section  4.3.1  recall  the  following  equation 

{M}  =  [h  ]{p}  (4.3.1.16) 

Jr 

and  from  Section  4.3,2  we  have 

{m'}  =  [h^]{p’}  (4.3.2.14) 
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Combining  (^1 . 3. 1. 16)  and  (^.3.2.14)  the  following  is  obtained, 

where  [H  ]  is  defined  as  [HP]  on  Page  15  of  Reference  ?• 

P 

{M*}  =  [H*]{p’‘}  (i|.3.3.1) 

In  the  same  manner  from  Equation  (^.3.1.12) 

W  =  LyJ[Q„]{M)  (i4.3.1.12) 

and  from  Equation  (^.3.2.10)  (the  axial  contribution) 

W  =  LBJ[q']{m'}  (i*.3.2.10) 

t 

Combination  of  the  [Qp]  and  [Qp]  matrices  yields  the 
matrix  [Qp]  which  is  defined  as  [QP]  on  Page  12  of 
Reference  7. 

The  final  pressure  load  vector  can  now  be  expressed  by  the 
following: 

{Pp}  =  [H]'^[Qp][Hp]{p*}  (i<.3.3.2) 

For  the  special  case  of  the  core  element  (Figure  2), 
the  pressure  load  vector  is  of  the  following  form: 

{F*}  =  [H]'^[Q*][H*]{p*}  (^+.3.3.3) 

4 . 4  Prestrain  Load  Vector 

The  prestrain  load  vector  is  constructed  assuming 
uniform  distribution  of  prestrain  across  the  element.  The 
prestrain  contribution  to  the  total  potential  energy  is: 

J  LeJ[E]{e.}dV  (4.4.1) 

V 


31 


where 


{£}  =  [D]{y} 

and  [D]  is  defined  by  Equation  (4.2,3) 

From  Equation  (2.7)  we  have, 

iy)  =  [H]{q} 

Therefore  Equation  (4.4.1)  becomes 

=  J  LYJ[D]'^CE]{e.}dV 

V 

We  also  know  that 


dV  =  2'iTrdrdz 


(4.4.2) 


(2.7) 


(4.4.3) 


=  LyJ  2Tr 


Let 


If 

r  z 


[D]'^rdrdz[E]{e^} 


(4.4.4) 


1 — 1 
to 

s 

-// 

[D]rdrdz 

(4.4.5) 

We  know  that. 

r  z 

0 

r 

0 

rz 

0 

0 

0 

1 

0 

1 

r 

z 

rz 

0 

0 

0 

0 

r[D]  = 

0 

0 

0 

0 

0 

0 

r 

r* 

(4.4.6) 

0 

0 

r 

0 

r 

0 

rz 
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And  from  our  previously  defined  notation: 


/  / 


drdz  =  I 


VJe  have  the  following 


1 10 

0 

Ill 

0 

0 

0 

0 

Iio 

loi 

Ill 

0 

0 

0 

0 

0 

0 

0 

0 

0 

Iio 

I20 

0 

1 10 

I20 

0 

Iio 

0 

Ill 

(4.4.7) 


Substituting  back  into  Equation  (4.4.4)  we  have: 


=  LyJ  C5]'^[E]{e^} 


(4.4.8) 


Recalling  that 


(y)’’  = 


(4.2.14) 


and  transforming  to  grid  point  displacement  coordinate  we 
have: 


♦e  =  LqJ  [H]'^C53'^[E]{e^} 


(4.4.9) 


Substituting  into  the  Lagrange  Equation  and  taking  the  first 
variation  with  respect  to  the  displacements,  we  obtain  the 
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prestrain  load  sector: 


{F  }  =  [.H]'^[D]'^CE]{e.}  (4.4.10) 

and  for  the  special  case  of  the  core  element,' we  obtain, 

{F^}  =  LH]'^[D][E]{e.}  '  (4.4.11) 


where 


(F  }  =  LFeI’i*  FgZj,  Fgr^,  FgZ^,  F^z^,  F^r^,  F^z^  J  ( 4 . 4 . 12 ) 


and 


=  L?i^»  ^  (4.4.13) 

4.5  Thermal  Load  Vector 

'  I 

The  thermal  load  vector  is  a  special  case  of  the, 
prestrain  load  vector.  The  temperature  distribution  function 
employed  for  the  trapezoidal  cross-section  ring  is  assumed  as 

I 

follows: 

J 

,  I 

T(r,  0,  z)  =  ki  +  k2r  +  ksz  +  k^rz  (4.5.1) 

or 


{T}  =  [g.]{k} 


■(4.5.2) 
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where  [g]  has  the  following  form. 


,1 

r  ' 

z 

rz 

1 

r 

z 

rz 

[g]  =  . 

1 

r 

i 

z 

1 

rz 

1 

r 

z 

rz 

The  prestrain  load  vector  contribution  to  the  total 
potential  energy  'may  be  written  as  follows:. 


The  initial  strain , vector  can  be  written  as, 
.  •  *  '  '  «  * 


{e.}  =  TLop,  ttg,  a.j,,  oj  -  {t5.} 


LeJCE]{e^}dV 

■  (4. 5., 4) 

V 

J 

From  our-  pz^evious  notation 

we  know  the  following:  . 

:  J 

{e}^  '  = 

; 

LyJOJ]’’ 

(4.2.4). 

'  (y)  ■  = 

t 

[H.]{g}'  ,  .. 

.(2.7) 

'  dV  = 

2Trrdrdz 

.  (^.5.5) 
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Upon  substitution  into  Equation  (^*.5.^)  obtain 


«!>£  =  2TrLqJ[H]Y  J 

r  z 


Rewrite  [E]{Ta}  as  [e“]{T} 


[D]'^[E]  {Ta}rdrdz 


where  the  coefficients  of  thermal  expansion  have  been 
multiplied  into  the  [E]  matrix. 


Equation  (4.5.6)  now  becomes. 


4'g.  =  27rLqJ[H]^J  J 
r  z 


rCD]'^CE^*]{T}drdz 


From  Equation  (4.5.2)  we  know  that 


{T}  =  [g]{k} 


♦g  =  27rLqJCH]'^J  J  r[D]'^[E“][g]drdz{k} 


r  z 


Define 


J  J rCDl^CE^lCgldrd 


z  as  [Q] 


r  z 


Then  Equation  (4.5.8)  becomes 


(4.5.6) 


(4.5.7) 


(4.5.2) 

(4.5.8) 


=  2iiLqJ[H]nQ]{k} 


(4.5.9) 
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The  {k}  vector  can  be  written  in  terms  of  the  grid  point 
temperatures  in  the  following  manner 


t 

I  L 


! 

k  i 


I 


I  I 


Ei 

! 

I  ^ 

%  P 


u 


r  I 

&  { 


I  ^ 


t 

r 


•  f 
I 


{<}  =  [h]{T‘} 


(4.5.10) 


where 


{k}^  =  Lki,  kz ,  ks,  k4j 


{TM  =  L(Ti-To),  (T2-T0),  (T3-T0),  (T4-To)J 


To  =  Equilibrium  Temperature  Of  Structure 


(4.5.11) 


(4.5.12) 


[h] 


-1  _ 


1 

ri 

Zl 

rizi 

± 

Tz 

Z2 

r2Z2 

1 

rz 

Z3 

r3Z3 

1 

1*4 

Z4 

r4Z4 

(4.5.13) 


and  [h]  has  the  following  form: 


[h]  =  i 


r2z;4(r3-r:,)  -r  1Z4 (r 3-ri, )  r4Zi(r2-rx)  -r3Zj(r2-ri) 

-Z4(r3-ri,)  Z4(r3-r4)  -Zi(r2-ri)  Zi(r2-ri) 
-r2(r3-r4)  ri(r3~r4)  -r4(r2-ri)  r3(r2-ri) 
(r3-r4)  -(r3-r4)  (r2-ri)  -(rj-ri) 


(4.5.14) 
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where 


A  =  (r2-ri)(r3-r.,)(z4-zi) 


Substitution  into  Equation  (^.5.9)  yields, 

=  2TrLqJCH]'^[Q][h]{T*}  (4.5.15) 


Substituting  into  the  Lagrange  Equation  and  taking  the  first 
variation  with  respect  to  the  displacements,  the  thermal  load 
vector  is  obtained: 


{F^}  =  27rCH]'^[Q][h]{T'} 


(4.5.16) 


For  the  special  case  of  the  core  element  the  thermal  load 
vector  is  of  the  following  form: 


{F^}  =  2TT[H]'^CQ][h]{T’} 


(4.5.17) 


i} .  6  Gravity  Load  And  Centrifugal  Force  Vectors 

The  work  done  by  the  acceleration  of  gravity  on  the 
Qisplacements  can  be  written  as: 


Work 


/ 


pGwdV 


(4.6.1) 
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where 


dV  =  2Trrdzc’r 

G  =  Acceleration  of  gravity 

p  =  Mass  density  of  the  material 

in  question 

w  =  Assumed  displacement  mode 

shape  in  the  z  direction. 


Work  =  2TrpG  J J'cBir  +  +  Barz  +  B4r^z)drdz 

r  z 


(ii.6.2) 


As  before  denote 


IF 

r  z 


z^^drdz  as  I_ 


pq 


(4.2.8) 


Work  =  l_Bi,  Ba,  Ba,  B^J  2TrpG 


110 

120 

111 

1 21 

V  / 


(4.6.3) 


Rewriting  the  work  equation  with  respect  to  all  the 
field  coordinate  degrees  of  freedom  we  obtain 
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V/ork  =  L*^!,  OLzy  as,  04,  61,  62,  03,  BuJ  2TrpG 


(J4.6.4) 


or 

Work  =  LyJ  {Fq}  (4.6.5) 

LYj  =  a2,  as,  a4,  0i,  02,  03,  04j 


~  »  I20  ,  In,  l2i_J 

(4.6.6) 

Remembering  that 

{y>  =  CH]{q} 

(2.7) 

lY>-  =  LqJ  chi'' 

C4.2.14) 

Work  =  LqJ 

(4.6.7) 

Upon  taking  the  first  variation  of  the  Work 
Equation  with  respect  to  {q}  ,  the  Gravity  Load  Vector 
{Fq}  is  obtained 

{Fq}  =  (4.6.8) 


4o 


The  Centrifugal  Force  Vector  Is  determined  as  follows; 
The  external  work  done  by  the  centrifugal  force  on  the  displacement 
can  be  written  as  follows: 


Work 


pw^rudV 


(4.6.9) 


where 


dV  =  2Trrdzdr 

ui  =  Natural  frequency  (rad/sec.) 

p  =  Mass  density 

u  =  Assumed  displacement  mode 

shape  In  the  r  direction. 

(4.6.10) 

(4.6.11) 


Work  =  27rpaj^  j  (u)r^dzdr 


Work  =  27rw*P  +  a2r*  +  asr^z  +  ai,r*z)dzdr 


Denote: 


J  J  r^z’^dzdr  as 
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where 

T  I 

{y}  =  L®i»  ®2»  p2»  P3»  ^'♦—1 

{Cq}*^  =  2Trp(u^Ll2 , 0 ,  l3,o,  l2,i,  l3,x,  0,  0,  0,  0,J  (M.6.15) 

Substituting  the  appropriate  transformations  into 
the  work  equation  and  taking  the  first  variation  with 
respect  to  {q}  we  obtain  the  Centrifugal  Force  Vector  {C^}. 

{Cq}  =  [H]'^{Cq}  (4.6.16) 


42 


4.7  Stress  Matrix 


The  element  stresses  are  given  by  the  following: 


ia)  =  [E]{e}  -  [E]{S} 


Recall  that 


(4.7.1) 


{£} 

— 

[D]{y} 

(4.2.2) 

r 

0 

1 

0 

"i 

0 

0 

0 

0 

{y} 

^0 

l/ri 

1 

r . 

0 

0 

0 

0 

{ 

>  - 

1 

0 

0 

0 

0 

0 

0 

1 

^i 

^  ^rzj 

0 

0 

1 

^i 

0 

1 

0 

(4.7.2) 

Prom  the  [D]  matrix  it  is  seen  that  the  stress  can  be 
evaluated  at  Nodes  (1)  through  (4)(i.e.  i  =  1-4) 
We  also  know  that 


{y}  =  CH]{q}  (2.7) 

Equation  (4.7.1)  can  be  written  in  the  following  manner: 

{o}  =  [E][D][H]{q}  -  [E]{a}  (4.7*3) 

And  for  the  core  element 

io)  =  CE][D][Hj{q}  -  rE]{a}  (4.7.4) 
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Denote 


CE][D][H] 

or  as  [S] 

CE][D][H] 


and  the  remaining  stress  contribution 


CE]{5}  as 


Equation  (^.7.3)  now  becomes 


{0}  =  [S]{q}  - 


where 


(4.7.5) 


(4.7.6) 


(4.7.7) 


(4.7.8) 
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V 


Kinetic  Energy 


It  is 
rotational 
energies . 
form: 


assumed  in  writing  the  element  kinetic  energy  that  the 
energies  are  small  compared  with  the  translation 
The  kinetic  energy  function  then  takes  the  following 


$2  =  m/2  J'iqMJ  [I]{q(m)}dV 


V 

The  assumed  displacement  modes  are  of  the  form 


LqJ  =  Lu,  wj 


(5.1) 


(5.2) 


Therefore 


or 


*2 


(5.3) 


(mil  u^  +  m22  w^  )dV 


(5.^) 


It  is  now  profitable  to  examine  the  form  of  the  assumed 
dispalcemsnt  modes  u  and  w 

u(r,0,z)  =  oi  +  a2r  +  aaz  +  a4rz 
w(r,e,z)  =  $1+02  +  63Z  +  3^rz 


(2.1) 

(2.2) 
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J 


Upon  examination  of  the  Kinetic  Energy  Function  (Equation  5.^) 
it  is  seen  that  all  terms  are  of  the  form: 

(a+br+cz+drz)  (e+fr+gz+hrz)  which  may  be  written  as 


La,  b,  c,  dj 


— 

— 

/  \ 

1 

r 

z 

rz 

e 

r 

r^ 

rz 

r^z 

z 

rz 

z" 

rz^ 

g 

rz 

r^z 

rz^ 

r^z^' 

h 

' 

^  / 

(5.5) 


For  the  element  of  volume  we  can  write 

dV  =  2Trrdrdz 


(il.2.6) 


Upon  multiplying  r  into  the  above  matrix  we  obtain  the  matrix 


which  is  to  be  integrated  over  the  volume 

r  r^  ,  rz 

r^z 

‘  » 

i 

r^  r*  r^z 

rz  r^z  rz^ 

3?  ^  Hi 

r^z^ 

.  (5.6) 

r^z  r®z  r^z^ 

r^z^ 

1  , 

Recalling  from  Equation  (4.2.8)  that 

, 

I  =  f  f  r^z'^drdz 

pq  Jj 

CO 

C\J 

.=r 

r  z 


The  kinetic  energy  function  can  now  be  assembled  in  matrix 
form  as : 


4)2  =  1/2  (_Yj  [M]{y) 


(5.7) 
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I 


Recalling  that 

iy)  =  CH]{q)  (2.7) 

(  ’  • 

and  \  . 

I 

■  '  iy)'  =  [H]{q)  ■'  ‘  C5.8) 

and' 

■  {y}*^  =  .q’[H]'^  ■  (5.9) 

Then 

I 

.  } 

$2  =  l/2LqJC'H]'^[M][H]{q}'  '  '  C5.10) 

Substitution  into  the  Lagrange.  Equation  and 
differentiating  once  with  •respecf  to  time  yields  the 
consistent  mass  matrix.  ■  ' 

•  CM]  =  [H]^[M][H]  '  (5.11) 

I  » 

For  the  special  case  of  the  core  element  the  consistent  mass 
matrix  is  given  as  follows: 

LM]*  =  [H]'^[M][H]  (5^2) 


U7 


•*  ^  -  ‘V  «.  •*  .  .-•  ->> 


LIST  OP  SYMBOLS 


CD] 

[E] 

{Fp} 

{F^} 

{Fq} 

CF^} 

CH] 

CH] 

[I] 

[K] 

CM] 

CS] 

T 

U 

W 

Lg] 

Ch] 

m 

P 

{q} 

r 

{/4) 


Centrifugal  Force  Vector 
Strain  Displacement  Coupling  Matrix 
Material  Properties  Matrix 
Pressure  Load  Vector 
Thermal  Load  Vector. 

Gravity  Load  Vector 
Prestrain  Load  Vector 
Transformation  Matrix 
Transformation  Matrix 
Identity  Matrix 
Stlffeness  Matrix 
Mass  Matrix 
Stress  Matrix 
Temperature 
Strain  Energy  Density 
Work 

Transformation  Matrix 
Transformation  Matrix 
Mass  Coefficient 
External  Pressure 

Displacement  vector  referenced  to  grid 
points 

Radius j  System  Coordinates 
Thermal  Stress  Vector 
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LIST  OF  SYMBOLS  (continued) 

u  Element  displacement,  r  direction 

w  Element  displacement,  z  direction 

z  Axial  Coordinate 

a  Field  Coordinate  Displacement  Degrees  of 

Freedom  Corresponding  to  displacement  in 

u  direction. 

a  Coefficient  of  Thermal  Expansion 

B  Field  Coordinate  Displacement  Degrees  of 

Freedom  Corresponding  to  displacement  in 

w  direction. 

e  Strain 

{y}  Vector  of  Combined  a  and  B  field 

coordinates 

K  Field  coordinate  Degrees  of  Freedom  for 

temperature  distribution  function 

V  Possion’s  Ratio 

a  Stress  Component 

u  Natural  Circular  Frequency  (rad/sec.) 

p  Mass  Density 

♦p  Potential  Energy  Function 

♦k  Kinetic  Energy  Function 
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C.  QUADRILATERAL  PLATE  ELEMENT 
I.  Introduction 

The  formulation  of  the  quadrilateral  plate  discrete  element 
described  is  derived  from  and  is  mathematically  consistent  with, 
the  formulation  described  in  Reference  12.  Tne  addition  of  this 
particular  element  serves  to  add  additional  capability  in  the 
analysis  of  shell  structures,  particularly  when  instability  analyses 
are  to  be  performend. 

A  detailed  derivation  ■*s  presented  for  the  force  displacement 
properties  of  an  orthotropic  quadrilateral  thin  plate  element 
exhibiting  membrane  and  banning  behavior.  Included  in  these 
relationships  are  terms  for  stiffness,  stress,  thermal  stress,  and 
incremental  stiffness. 

For  the  quadrilateral  plate  element,  orthotropic  material 
mechanical  properties  are  defined  by  four  parameters:  E  ,  E 
and  Gxy  where  ^^xy  Poisson's  ratio  of  the 

contraction  in  the  y  direction  to  extension  in  the  x  direction 
due  to  a  tensile  stress  in  the  x  direction.  There  is  another 
Poisson's  ratio  -^yx  similarly  defined,  which  is  related  to  the 
other  material  properties  through  the  identify  E^^^x  “  ®y'^y* 
Since,  in  general,  none  of  the  sides  of  an  element  will  correspond 
to  a  principal  axis  of  orthotropy,  all  relationships  are  derived 
for  an  arbitrary  orientation  of  the  element  in  an  x-y  plane  in  which 
the  X  and  y  axes  are  parallel  to  the  principal  orthogonal  directions 
of  the  material. 

Techniques  for  deriving  the  desired  force-displacement 
relationships  are  described.  The  derivations  make  use  of  the  matrix 
statement  of  Castigliano's  Theorem.  As  shown  in  Reference  13,  this 
leads  to  the  stiffness  matrix  being  defined  as  the  product  of  three 
matrices,  each  containing  simple  terms. 
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Development  of  Linear  Elastic  Membrane  Stiffness  Matrix 


/ 


I. 


For  the  quadrilateral  plate  element, 
the  stresses  are  assumed  to  be  given  by 


shown  in  Figure  III -3, 


X 

'xy 


=  a^^  +  a^y 


=  ar 


(III-l) 


This  assumed  stress  pattern  was  adopted  in  Reference  13  by  Turner, 
et  al,  in  the  derivation  of  the  stiffness  matrix  for  an  isotropic 
rectangular  plate.  Equations  (III-l)  satisfy  the  equilibrium 
requirements  of  a  differential  element  and  can  be  operated  upon 
to  yield  compatible  displacements.  The  complete  element,  in 
attachment  to  other  elements  of  the  system,  does  not  satisfy 
equilibrium  and  compatibility  at  all  points  along  the  juncture 
line  however.  Evidence  (Reference  5)  from  plane  stress  analyses  has 
shown  that  the  consequence  of  these  shortcomings  is  a  stiffer 
idealization,  but  not  as  stiff  as  would’ be  obtained  with  use  of 
the  linear  edge  displacement  assumptions  proposed  in  Reference  6. 

The  matrix  statement  of  Castigliano • s  Theorem,  applicable 
to  the  derivation  of  the  quadrilateral  element  force-displacement 
properties  is  written  as 

[k]  =  [b’^J  [cJ  [b]  (III-2) 

where  [t]  is  the  desired  matrix  of  element  stiffness  coefficients. 
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P-TPiURE  III- 3  GFO'IETRY  OP  ARBITRARY  QUADRILATERAL  PLATE  ELEMENT 


i 


'i 


To  obtain  the  matrix  [5]  j  we  first  establish  the  strains 
for  an  orthotropic  material  in  accordance  with  Hooke's  Law  and 
considering  the  possibility  of  initial  strains  as  follows: 


£'x  ~^xy  ^  "^^x 

^y  “  ET"  ~^yx^  ^  ‘^^y 


(III-3) 


y  =  ^ 

®xy  5]^  *xy 


i  i 

where  *  ^y 


and  are  the  initial  strains. 


The  constants,  a,  in  the  assumed  stress  pattern  are 
introduced  into  the  strain  expressions  by  substituting 
Equations  (III-l)  into  the  above  equations  to  yield 


^x  “ 

^  [ 

(  ®l'^y  ^3^ 

--^xy  V  +  “gyj 

[ 

V 

^  ®3"^yx  ^1^ 

-Ax  *2''  ^  v] 

Utilizing  the  basic  relationships  that 

g  ,  it  is  possible  to  determine  the  linear  displacements 

^  y  y 

of  the  quadrilateral  by  integrating  the  strain  expressions 
^Eq.  IJ.I-4)  with  respect  to  the  appropriate  variable.  Thus, 

l€^  dx 


J  X 


C-  [(«i-^xy“3)  y  '  +  «2y=' +  fi(y)]  + 


(III-5) 


=Je/y 


2 


The  functions  of  integration,  f2^(y)  and  f2(x),  can  be  evaluated 


from  the  shear  strain  definition 


/  = 

"xy  d  y 


dx 


Substituting  Eqs.  III-4  and  -5  into  this  shear  strain  definition 
gives 


»5  .  y  ^ 

^  'xy  -  -Eg- 


1  4  fi(y) 


1  d  f2(x) 


ET"  ET  -aST 


(m-6) 


To  separate  the  variables,  Eq.  III-6  is  rearranged  as  follows 


a  f,(y)  ,  VI  .S  ./  ^l(iii.6a) 

^x  dy  ®y  L  ^  ^x  ®xy  -* 

The  functions  f|^(y)  and  determined  by  letting  each 

side  of  Eq.  III-6a  equal  the  constant  ®6  and  integrating  to 
yield  x 


^l(y)  =  "  IT  *6^ 


(III-7) 


^2  i-  V  "  +  ®y  ^xy  * 
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Substituting  Eqs.  III-7  into  Eqs.  III-5  and  evaluating  the 
resulting  expressions  at  the  corner  points  1,  2,  3  and  4  of 
Figure  1  produces  the  relationship  for  the  corner  displacements 
as  shown  as  Eq.  III-8  in  Pig.  Ill- 4.  The  square  matrix,  including 

the  coefficient  1 _  ,  on  the  right  side  of  Eq.  III-8  is  the^Bj  matrix 

' 

To  develop  the  (|c]|  matrix  we  first  need  the  expression  for 
strain  energy  for  orthotropic  plane  stress.  In  terms  of  stress, 

- 2GL-  cr  <r  dx  dy  (III-9) 

E  ^  y  ^xy 


this  can  be  written  as 


Substituting  the  assumed  stress  l*Jinctions,  Eqs.  III-l,  into  the 
energy  equation  and  expanding  produces 


2  2 


2  2 


W  ^  ^  ^  ^  tm. 

(a^  +  2a^  agy  +  ag  y  )  +  ^  (^3  +  Sa^  a^^x  a^^  x  ) 


a. 


(  a^^a^  +  a^a^^x  +  ag  a^y  +  a^a^j^  xy)  +  j  dA  (III-IO) 


] 


The  derivatives  of  the  strain  energy  with  respect  to  the  constants 
a^^,  ag,  . a0  are: 


E. 


<)ai 

11 

I^Aai 

J  U 

^ag 

h 

“  IT 

X 

[Vi 

2 

J  u 

aas 

h 

"  IT 

X 

1 _ t 

(■’V'' 

Ju 

^  a4 

h 

"  ir~ 

X 

i _ 1 

E. 


x®4j 


] 


'•xy“4 


fIlI-11) 


‘-^yx^xy®2 


+  +  »lt)] 


^t?^y^?J'j^ff'fflfe;aTJ*y;yiBTWa'Wpyt'ap>v»iew('»>i  ywj  ■»  <■  ^.»»^^  v» 


d  a^ 

_  h  fv 

•>] 

du  _ 

^  a„ 

7 

<)u 

Jag 

=  0 

where 

A  = 

/  dA 

and 

i  J 

I  X  y 

=  J 

^  i  J 
X  y 

A 

dA 

Prom  the  above  we 

have 

•« 

A 

■^x/ 

2 

2 

^xy  ly 

^y^xy 

^  U 
^ao 

3 

■'V' 

- 

^y^y 

E 

3  A 
IT  " 
y 

E 

IT  ^x 
y 

‘  h 

"■^y^x 

^y^xy 

E 

^  T 

^Ix 

®x  j 

V  " 

du 

0 

0 

0 

0 

du 

0 

0 

0 

0 

du 

0 

0 

0 

0 

^  a™ 

7 

0 

0 

0 

0 

(III-ll  Cont.) 


V 


0  0  0  I  fa 

0  0  0  I  I  a, 

0  0  0 

0  0  0  I  la, 

0  0  0 

0  0  0 

0  0  0  1  I  a, 

0  0  0  I  I  a< 

(III-12) 

The  square,  symmetric  matrix,  including  the  coefficient  |i-  , 
on  the  right  is  the  Matrix  fcj  .  * 

Utilizing  Matrix  [bJ  from  Eq.  III-8  and  Matrix  [cjfrom 

Eq.  III-12,  the  element  stiffness  matrix,  Lk]  ,  is  obtained  from 
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Eq.  III-2. 


8 


*  ^  •  »'*  4,  y. 


III.  Geometric  Properties 

In  the  development  of  Matrix  [c]  in  the  previous  section, 
a  group  of  I„  J  terms  resulted  from  integrating  the  energy 

X  j 

expression.  These  terms  are  of  the  form 


=jJxV  dy  dx  =J^  xV*^  dA 


(III-13) 


and  therefore  are  analogous  to  area  moment  relationships. 

Tvi  *1  o  1  4-1  11*  o4"*l  nrr  4-Via  T  ^  ^  4*Aiimo  4 


In  explicitly  formulating  the  I 


terms  it  was  con¬ 


venient  to  use  previously  determined  ^  terms  with  reference 
to  the  ^  coordinates  (see  Pig.  III-3)  and  transform  these 
to  the  x-y  coordinate  system.  The  coordinate  trems formation  from  . 
the  x-y  system  into  the  coordinates  for  any  given  point  is 
expressed  by: 


cos  6  sin  Q 
sin  6  cos  Q 


m 


(III-14) 


where  0  is  the  angle,  betiieen  the  x  and  ^  axes  and  is  defined  by 


e  =  tan 


(III-l4a) 


Alternatively  the  coordinate  transformation  from  the  system 

to  x-y  system  is  given  by 

fxj  ^Fcose  -slne]w|  (III-l4a) 

/  yj  sin  0  cos  ©j/’J 

As  an  example  of  the  determination  of  's  in 

i  -t  *  ^ 

terms  of  I  f  ^’s  consider  the  first  moment  of  the  area  about  the 
y  axis,  I  .  For  this  moment 

X 

=  fa  X  dA 


i  i  . 


s  in 


=  j  (7^  cos  0  -  5sin  0  )  dA 

=  cos  Q  jjfdPi  -  sin  0j'^  dA 


(III-15) 


cos  0  - 


sin  0 
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All  the  other  I  ^  ^'s  are  determined  by  a  similar  procedure 
X  y 

and  are  stated  as  follows: 

cos^G  -  cos  0  sin  6  +  sin^G 

ly  ~  ^7^  G  +  cos  0 

ly^  =  1^  sin^G  +  21;^^  sinG  cos  G  +  cos^G 

Ixy  =  -  I§^)  sin  0  cos  0  +  (cos^G  -  sin^G) 


(III-15a) 


The  area  and  moment  properties  about  the  7^  and  ^ 
axes  were  determined  by  direct  integration  within  the  proper  limits 
to  yield  the  following  expressions 

A  =J'dA  =  ^  (^^4^4  +  (^3“  ^4)  (f3+^!4)  “ 

“  5:  f  (^4^3  “C3/?4)  (^3+^14)  +^3^2 

^  ll  [  %  +  ^3^4'^^4^)  ■‘■^3^2  (^3^  '^^3^2  ‘‘'^2^0 

If  +  (?3“^4)  (^3^+ ^3^4  +^4^)  "  (^^2)^3^]  (III-16) 

If^  “  X2  |^^4^^4  (^3~^4^  (^3  **■  ^4  )  (^3'*^4^  “  (^3"^2^ 

I,  _1  [3^21,2  3  (^3-^4)  (^.3^  -^4^  ) 

^71^-  ^  L3f4  ^4  -  - - 

^  3  (^3“^4^  ^?^3  '^?3^4  "^^4  ^  ^^4^3  "^3^4^ 

(-^3  “^4  ) 


IV.  Development  of  Initial  Force  Terms 

The  Initial  force  terms  are  derived  from  a  consideration 
of  the  element  corner  forces.  Prom  Castigliano's  Theorem,  the 
corner  forces  are  expressed  by 


W- 


(III-17) 


Noting  that  the  matrix  is  the  set  of  influence  coeffi¬ 
cients  in  the  direct  relationship  between  the  constants  a  and 
the  displacements  /(=  u,  v),  it  follov/s  that 


(III-18) 


Also  the  vector  already  been  determined  in  Equ.  III-12, 

Sect.  III-A,  as  =  [c]  [a|  .  Thus,  we  can  rewrite 

Eq.  III-17  as 


I  p]  =  [  B-l]  [cj  {aj 


(III-19) 


The  initial  strains  are  introduced  into  the  force 
expression  by  first  solving  Eq.  III-8  for  [a}  thusly 


W-  [•]■'  M  ■  W-  ftl; 


(Ill-8a) 


and  substituting  this  relationship  into  Eq.  III-19  to  give 

f.j .  [.-•]  ’  [o]  ’  [0]  [.]  ■'  (I;!; . , 


XV  ^ 

(III-20) 


The  first  term  on  the  right  hand  side  of  Eq.  II 1-20  represents 
the  corner  forces  due  to  displacement,  i.e.,  the  forces  required 
to  Induce  the  deformations  u  and  v  elastically.  The  second  term 


60 


(III-21) 


yields  the  initial  forces,  ^  P^J.  Thus 

^[s-fCc]  [bJ  -1  ^  ^  j  (III-21) 

T 

Since  from  Eq.  III-2  [k]  =  [b'^J  [c]  [bJ  the  initial  forces 

are  simply  the  product  of  the  usual  element  stiffness  matrix  and 
the  column  of  node  point  initial  displacements.  Hence 


X  *3 


[kJ 


(m-22) 


^■yVa  +  ^xy^*3 

€  V  t  r  \ 

y  xy 


For  the  case  of  a  change  in  the  termperature  of  the 
quadrilateral  element  equal  to  T,  the  initial  strains  are  given  by 

^  i  _  c»<  T 

X  X 

f  1-  T  (111-23) 

y  y 

^xy"=  ° 

where  <<  and  «»<  are  the  coefficients  of  thermal  expansion 
X  y 

in  the  X  and  y  directions  respectively. 


V.  Stress  Equations 

The  stress  equations  are  derived  by  first  evaluating  the 
assumed  stress  functions,  Eqs.  III-l,  at  corners  1,  2,  3  and  4. 
This  procedure  yields  the  following  relationship  for  the  corner 
stresses: 


<^x 


■  \ 
^1 

^6 
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4 
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‘■xyg 

'T 

^xy3 

<^Y4 
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0 

0 

0 

0 

0 

0 

0 

0 

0 
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0 

0 

0 

0 

0 

0 
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0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 
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In  a  more  concise  manner  Eq.  III-24  may  be  written  as 

[•]  H 


“3 

®4 

a^ 


7 


(III-24) 


(III-24a) 


By  substituting  Eq.  III-8a  from  Sect.  III-C  into  Eq.  III-24a, 
the  stresses  may  be  written  in  terms  of  the  displacements  and  initial 
strains.  This  substitution  yields 


51  =  w  [b]'^  W"' 


C  ^  y  +  r  ^ 

»  •  ir 


xy 


(111-25) 


62 


The  first  term  on  the  right  hand  side  of  Eq.  III-25  corresponds 
to  the  displacement  stresses^^J  ,  v*h'*le  the  second  term  represents 
the  initial  stresses, [<r  i}  . 

Letting  [s^y]  =[dJ  [b]"^  ,  the  displacement  stresses  are 


given  by 

/V„  A 

"I 

i 


u 


X3 

X4 


'L 


xy. 


't  ^ 

xy2 


‘C 


xy- 


'L  i 

xy4 

cr  i 

yi 

CT  ^ 
^2 


^1 

^2 

^3 

-^4 


/ 


^1 

^2 

^3 

.^4 


(III-26) 


-  ^ 


VI.  Development  of  Linear  Elastic  Stiffness  Matrix  (Bending) 

The  orientation  of  this  element  in  the  X-Y  plane  and  the  corner 
forces  and  moments  are  shown  in  Pig. 111,5,  As  before,  the  principal 
directions  of  orthotropy  must  be  parallel  to  the  X  and  Y  axes 

In  deriving  the  linear  elastic  stiffness  matrix  the  displace¬ 
ments  in  the  direction  normal  to  the  plate  midplane  are  assumed  to 
be  representable  as 


w  =  a^x-'  +  a^x  +  a-x  + 


^x  +  a^^y-^  +  a^y  +  a^y  +  a^-'y  +  agX  y 


+  a„xy  +  +  a^j^xy^  + 


(V-1) 


Where  a^, . a^^g  are  constants.  Since  twelve  force-displacement 

equations  are  to  be  derived  (three  degrees  of  freedom  -  0  ,  0  ,  w  - 

X  y 

at  each  corner  point),  twelve  independent  parameters  appear  in  the 
assumed  displacement  function.  Secondly,  all  possible  single  terms  or 
products  to  the  third  degree  are  included;  the  resulting  polynomial 

O 

is  "geometrically  symmetric",  e.g.,  corresponding  to  the  agX  'y  term 

p 

there  is  the  terra.  Finally,  Equation  (1)  satisfies  the 

differential  equation  of  equilibrium. 

To  obtain  the  desired  stiffness  matrix,  the  following  matrix 
product  is  effected; 


([b]  [oj  [b] =  [kj.] 


(V-2) 


To  formulate  the  matrix  bJ  ,  we  first  establish  the  slopes  0  and  0 

X  y 


from  the  deflection  function  (Equation  V-1)  by 

0^  =  =  Sy^a^^  +  2ya^  +  + 

^  y 


®y  = 


2  ^2 

=  3y  a^i^  +  2ya^  +  ^  +  *  ®8  ^*9 


(V-3) 


+  3xy  h  Sxyaj^j^ 

2  2 

=  -3x  a^^  -  2xa2  -  *3  -  3x  ya,^  -  2xyag  -  ya^ 

3  2 

-  y  am  -  y  a-n 


(V-4) 


The  slopes  and  the  deflection,  w,  at  the  four  corners 
of  the  quadrilateral  are  evaluated  to  yield  the  relationships  shown  as 
Equation  (V-5)  in  Figure  (V-  2  ) .  The  square  matrix  on  the  right 
hand  side  of  Equation  V  5  is  the  [ B]matrix. 

The  [  j  matrix  is  determined  by  first  considering  the 
flexural  strain  energy.  In  terms  of  the  out-of-plane  deflection 
w,  the  strain  energy  is  given  by 


where: 


^  ^  '  '^xy  '^yx 


is  shown  in 


Figure  III-7. 


(V-7) 
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PIOUREIII-7  LC.J  MATRIX  FOR  QUADRILATERAL  FLEXURAL  ELQ(E::T 


By  calculating  the  inidicated  partial  derivatives  of  v  from 
Equation  V-1,  and  substituting  into  Eq.  V-6,  the  following  expression 
for  the  energy  in  terms  of  the  constants,  a,  is  obtained 


^  1/  +  24xa3^a2  -}■  72x^ya3^a^  +  24iya^aQ 

^  P  P  P 

+4a2  +  24xya2a^  +  8y®2^3  y‘»^' 

+24x  y^agag  +  4y^ag^) 

E 

+  ^  f  36y^8iu^  +  24ya2^ac  +  72xy^a2^a^Q 

X  ' 

+  24xya2^a^j^  +  4a^^  +  24xya^aj^Q  +  Sxa^^a^^j^ 

+  36xVa3^Q^  +  24x^ya3^Qaj^3^  +  4x^a^3^^j 

+  2/^y  ^  (  36xyaj^a2^  +  12ya2a2^  +  36xy^a2^ay 

+  12y^a2^ag  +  12xa2^a^  +  4a2a^  +  12xya^ay 

2  2  2 
+  4ya^ag  +  SSx^ya^^a^^Q  +  12xya2aj^Q  +  36x  y  a^a^^^ 

+  ISxy^aga^^Q  +  ISx^a^^a^^^^  +  ^xaga^j^  +  12x^ya^aj^j^ 

+  4xyagaj^^) 

4Gvv  M  /  4  2  r-  2 

+  — M. -  ^  9x  By  +  12x-'ayag  +  6x  a^a^ 


4GvV  ^ 


+  12x^ayag  +  Sx^a^a^ 


+  iSx^y^a^aj^Q  +  12x^yayaj^j^  +  4x^aQ^  +  4xagag 
+  12xy^  aga^Q  +  8xyaga^^  +  a^^  +  6y^a^a^Q  + 


Next,  the  strain  energy  is  differentiated  with  respect  to  each  of 

the  constants,  a^^, . a^2‘  example,  the  derivative  with 

respect  to  a,  is  _ 

j.  f  ' 

d  /  >.  2  ^2  N 

=  -jIjj—  LC36x‘^aj^  +  12xa2  +  36x  ya^  +  12xy  agj 

+  2y<^xy  —  I  iSxya^^  +  6xa^  +  l8x  ya^^ 


+  6x^a, 


dx  dy 


(V-9) 
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Letting  JJ  x  dx  dy  =  I^,  or,  in  the  general  form  letting 
dx  dy  =  »  Eq.  V-9  is  rewritten  as 

■3 

Up  Ev  ^  r  p  p 

7^  =  -W  L  36  Ix^l  +  36  l/ya^  +  12  I^yag 


t^xy  if  lxy®4  Ix®5  +  36  ya^^Q 

+  12 


(V-IO) 


The  development  of  all  terms  will  he  presented  in  detail 

X  y 

in  Section  C. 

Three  of  the  derivatives  are  zero;  namely. 


=  0 


In  matrix  notation  the  twelve  partial  derivatives  of  the  energy 
with  respect  to  the  constants  may  he  expressed  as 


-h]  • 


The  square,  symmetric  matrix  C^.  J  is  given  in  Figure  III-7. 

The  |^B]and  matrices  can  now  he  used  to  determine  the 
flexural  stiffness  matrix,  [  ,  for  the  quadrilateral  plate 

element  from  Eq.  V-2.  Clearly,  the  operation  fs-l]  '^is 

too  complicated  to  allow  for  an  explicit  formulation  of  the  [  K^J 
matrix.  It  is  Intended  that  the  formulations  of  the[B3  andfc^j 
matrices  he  stored  in  the  computer  for  use  each  time  a  stiffness 
matrix  is  to  he  evaluated.  Unfortunately,  it  has  been  found  that 
for  certain  geometric  proportions  the^B^  matrix  will  he  singular. 
A  means  for  predicting  this  singularity  and  circumventing  it  is 
discussed  in  Section  X. 
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VII. 


Development  of  the  Incremental  Stiffness  Matrix 

The  incremental  stiffness  matrix  f n^  is  derived  through 
application  of  the  matrix  triple  product 

■^)  ^  [“n]  [b]  =  [o] 


The 


[b] 


(V-12) 


matrix  is  the  same  as  that  which  has  oeen  discussed 


above.  The  [  C^j  matrix  is  defined  in  a  similar  manner  as  before, 
but  not  the  appropriate  "energy"  integral  is 

“n  =  [  “x  (7^)  "  +  "ydr)  "  ^  (jl)  (jj) 

In  developing  the  ^  matrix,  it  is  assumed  that  the  inplane 

forces  N„,  N  and  N  are  constant  throughout  the  element.  Since 
X  y  xy 

the  normal  inplane  stresses,  o'  and  o'.,  are  actually  not  constant 

X  y 

(see  Chapter  III),  the  inplane  forces  are  taken  as  the  average  of 
the  edge  forces  occurring  at  the  four  corner  points  of  the 
quadrilateral  plate  element.  Thus 


^<^x,  '^‘^1 


^^1  ^2  ^^3  ^^4 


-} 

-) 


h 


The  shear  stress  is  constant  throughout  the  element  so  that 


\y  ^xy  ^ 


xy 


It  is  convenient  to  divide  the  energy  expression.  Equation 
V-13,  into  its  three  components  as  follows: 


U„  =  U 
n  n. 


+  U. 


n. 


+  U, 


n. 


1 


N  {{ ^ 
\j\  [TxJ 


xy 


dx  dy  +  ^ 


N, 


(iff)  "xyljl^ 


dx  dy 


(V-13a) 


72 


We  now  consider  each  energy  component  separately  and  obtain  a 
matrix  corresponding  to  each  component  (J  =  x,  y  or  xy) . 


(s) 


By  differentiating  the  assumed  deflection  function  (Equation  V-1) 
and  substituting  the  derivatives  into  Eq.  V-13a,  we  obtain  an 
expression  for  each  energy  component  as  follows: 

Un  =  ^  l^x^a^^ag  +  Sx^a^a^  +  iSx^a^^a^  +  lax^a^^ag 

P  p  O  P  P  P  P 

+  6x  ya^a^  +  6x  y^a^^^lO  ^  ®1®11  + 

+  12x'^ya2a^  +  8x  yagag  +  4xya2a^  +  ^xy^aga^^  +  4xy  a2®^li 

+  a^  +  6x  ya^a^  +  4xya2ag  +  2YB.^ei.^  +  2y-’a2a^Q  +  2y 

+  9x^y^a^^  +  12x^y^a^ag  +  Sx^y^a^a^  +  6x^y^a,^a^Q  +  6x^y^a,^aj^j^ 

+  4x^y^ag^  +  4xy^agag  +  W^aga^^  +  4xy^aQa^^  +  y^a^^ 

+  2y^a^a^Q  +  Sy^a^a^^^^  +  y^a^^^  +  +  y\l^] 

(V-l4) 


U. 


”y 


1 


Ny J  [9y^a4^  +  ISy^ai^a^  +  Sy^a^^a^  +  6x^y\a^  +  6x^y%ag 
+  exy^a^^a^  +  iSxy^a^a^Q  +  12xy3a2^a^j^  +  4y^a^^  +  4ya^ag 
+  4x^ya^a^  +  4x^ya^ag  +  4xya^ag  +  12xy^a^a^Q  +  Sxy^a^a^^^^ 

+  ag2  +  2x\a^  +  ax^agag  +  2xagag  +  exy^a^a^^  +  4xyaga^^ 


+  x^a^^  +  2x^a^ag  +  2x‘*a^a^  +  6x‘*y‘^ayaj^Q  +  4x‘*ya^aj^j^  +  x'^ag 


4„2. 


.4^  2 


+  2x'^aga^  +  Sx^y^aga^Q  +  4x^yagaj^j^  +  x'^a^'^  +  6x‘^y‘^aga 

+  4x^yaga^3^  +  9x^y^aj^Q^  +  12x^y3aj^Qa^j^  +  J  dx  dy 

(V-15) 


4x^: 


2.  2 


.2.2. 


10 
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^xylf  +  Sy^a^a^^  +  Sx^y^a^^a^ 

o  o  c  2l  P 

+  6xy^a2^ag  +  3y^a2^a^  +  Sy^a^^a^^Q  +  3y  a^^a^^j^  +  ^x  ya^^a^ 

+  4xya2a^  +  Sya^a^  +  6x^y~a^a^  +  4xy^a^ag  +  2y^aga^ 

4  ■R  2 

+  2y  a^a^Q  +  2y-’aga^j^  +  3x  a^ag  +  Sxaga^  +  a^ag 
+  SxSaga^  +  Sxyagag  +  yagBg  +  y\aj^g  + 

c  4  1^  2  k 

+  3^  SL^&r^  +  2x  aga^  +  x-^a^a^  +  3x'^ya^  +  3x  a^^ag 

3  2  4  3  2  3 

+  2x“'a2ag  +  x  a^ag  +  5x  ya^ag  +  2x"’yag  +  3x'^a^ag 

+  2x^a2ag  +  xa^a^  +  4x^ya^a^  +  3x^yaga^  +  xya^^ 

+  9x^y^aj^a^Q  +  6x^y^a2aj^Q  +  3xy^a2a^Q  +  lOx^y^a^a^^ 


Next  the 


+  Tx^y^aga^Q  +  4xy3aga^Q  +  3xy5a^Q^  +  Sx^ya^a^^^ 

+  ^x^ya2a^^  +  2xya2aj|^^  +  7x^y^a^a^^  +  5x^y^agaj^^ 

+  3xy^a^aj^j^  +  Sxy^a^^^a^^^  +  2xy^a^^^J  dx  dy  9(V-l6) 

5  jC  I  matrices  are  developed  by  taking  the  derivatives 

L  Jj 


of  the  respective  energy  component  with  respect  to  the  constants, 
a.  For  example  in  deriving  thej^C^  J  matrix,  the  differentiation 

of  Ujj  with  respect  to  a^^,  is 


*x  = 


<ra][" 


\  11  iSx^aj^  +  12  x\2  +  iSx^ya^  +  12x^yag 

2  p  2  3  2  2  T 

+  6x  ya^  +  6x  y  a^^^  +  6x  y  a^^^^  J  dx  dy 
\  f  ^^x  ®1  '*’  ^^x  ®2  ^^x  ®3  ^^x  y®7^^x  y^8 


^^x  y®9  **■  ^^x  y  ^10  ^^x  y  ®11  J 


(V-17) 
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As  noted  previously,  the  terras  are  developed  In  Section  C. 

jT 


All  the  other  derivatives  are  calculated  sirailarly  and 
lead  to  the  following  expressions 


The  matrices  jc  1  #  [c  1  ,  |  are  shown  in  Pigs.  III-8,  -9 

L  "xj  L  L  ^yJ 

and  -10  respectively. 

Since  jc  1=10  +  jc  |+  C  ,  the .  incremental 

L  "J  L  "xj  L  "yj  I  "xyj 

stiffness  matrix  may  be  stated  as 


(V-12a) 
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MATRIX  FOR  QHJADRIUVTERAL  FLEXURAL  ELWENT 


VIII.  Geometric  Properties 


li 


H 


The  derivation  of  the  I. 


terms  in  the 


[oj  an(j|c^^ 


matrices  are  presented  in  this  section.  These  terms  depend  only  on 
the  geometry  of  the  element  and  are  defined  as 


'xV 


y^*  dx  dy 


Since  the  terms  I^  where  is  an  axis  coincident 
with  one  side  of  the  quadrilateral  (See  Figure  III-3)  had  previously 
been  obtained  by  direct  integration,  it  was  convenient  to  express 
the  moments,  j  in  terms  of  the  already  known  's. 

The  transformation  to  the  coordinates  from  the  x-y 

coordinates  is: 


cos  6  sin  0 

-sin  0  cos  0 


(V-19) 


and  alternately 


cos  0  -sin  0 
sin  0  cos  0 


(V-20) 


where: 


sin  0  = 


cos  0  = 


To  illustrate  the  derivation  of  the  I  ^  ^'s  consider 

X  y 

p 

the  second  moment  of  area  about  the  y  axis,  I  .  By  definition 
=  J  (^cos  0  -  §  sin  0)^  dA 

=  f  (h^cos^0  -  2  ^^sin  0  cos  0  +  ^  ^  sin^  0)  dA 

t  ^  f^r  m 


cos^  0-21^^  Sin  0  cos  0  +  sin^  9 


(V-21) 


'avwiw*-; 


I 

K, 

I 

I; 

? 

r 

I- 


i 


This  procedure  was  used  to  determine  all  the  terms 

«y 

of  the  1^5*^ 's.  The  have  been  derived  elsewhere. 


The  area.  A,  of  the  quadrilateral  can  be  directly 
determined  by  adding  and  subtracting  triangles  and  trapezoids. 
In  both  coordinate  systems  the  area  is  given  by 


A  =  I  (  x^y^  +  x^y^^  -  xj^y^  -  x^yg  ) 

^  =  I  (^243  -^4  f3  ) 


(V-22) 
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IX.  Development  of  Corner  Thermal  Moments 

Although  the  corner  thermal  moments  can  be  formulated 
through  the  use  of  Castigllano ‘s  First  Theorem  (a  procedure 
consistent  with  the  derivation  of  the^K^j  and  [nj  matrices),  it 
is  considerably  simpler,  and  probably  sufficiently  accurate,  to 
obtain  the  thermal  moments  by  pro-rating  the  distributed  edge 
moments  to  the  corners.  This  direct  approach  was  used  here  and 
is  described  in  the  following. 


For  the  case  of  a  temperature  variation  through  the  thick¬ 
ness  of  the  plate  element,  the  thermal  moments  per  unit  length 


at  any  point  i  is  given  by 


M. 


M.r  = 


(V-23) 


where  §  is  a  thickness  coordinate  measured  positively  in  the 
positive  z  direction  from  the  neutral  axis  of  the  cross-section. 

In  deriving  the  corner  thermal  moments,  it  is  assumed 
that  the  distributed  thermal  moments  are  constant  throughout  the 
element  and  equal  to  the  arithmetic  average  of  the  moments  at 
the  four  corners  of  the  quadrilateral.  Thus 


V  ( v’  v'  ^  "  'vr*  > 


(V-24) 
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These  average  moments  are  distributed  around  the  edges  of  the 
quadrilateral  element  as  shown  in  the  following  sketch; 


y 


The  distributed  thermal  moments  are  concentrated  ("lumped") 
at  the  corners  of  the  element  by  assigning  one-half  of  the  total 
moment  along  an  edge  to  each  corner  bounding  the  edge.  For 
example,  the  corner  thermal  moments  at  Corner  1  are  obtained  by 


I  My  (yg  -  y^) 


(V-25) 
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The  corner  thermal  forces  Fj^  are  zero  so  that  the  thermal  moments 

z 

and  forces  are  expressed  in  matrix  notation  as 


By  differentiating  the  assumed  displacement  function  (Eq.  V-1), 
the  following  derivatives  are  obtained 


6  ^w 

dx^y 


=  6aj^x  +  2b.^  +  6a,^y  +  2a0y 


=  6a||^y  +  2a^  +  ^®^ll* 

=  +  2aQX  +  a^  +  Sa^^^y^  +  2a^^y 


(V-28) 


Substituting  Eqs.  V-28  into  Eq.  V-27,  and  evaluating 
the  expressions  at  the  centroldal  coordinates  (  x,  y)  yields 
the  following  matrix  equation  for  the  moments  and  forces: 


=  [-1  r 

I  • 


(V-29) 


The  ^gj  matrix  is  shown  in  Fig.  m-n. 

The  column  of  assumed  constants,  E*!*  V-29  is 

eliminated  by  solving  for  fa^^lin  Eq.  V-5,  giving 


:  =  [bJ 


(V-5a) 
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I 


i 

7 

X 


1^ 


are  the  average  distributed  thermal 


In  Eq.  V-32,  and 

moments  defined  by  Equation  V-24. 
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X.  Criterion  for  the  Singularity  of  the  Matrix 

Under  certain  geometric  conditions  the  [b]  matrix  may  he 
singular.  When  this  is  the  case  there  is  no  alternative  but  to 
revise  the  analytical  model  so  as  to  define  an  element  of  different 
geometric  proportions.  It  is,  of  course,  undesirable  to  permit 
a  complete  analysis  before  the  singularity  is  recognized.  Hence, 
the  following  criterion  for  assessing  singularity  or  111-condi¬ 
tioning  has  been  developed;  it  can  be  applied  by  means  of  hand 
computation  before  the  analysis  is  performed  or  Incorporated 
as  a  check  in  the  routine  for  the  computation  of  the  element 
force-displacement  relationships . 

The  derived  criterion  is  actually  the  result  of  an  attempt 
to  develop  an  explicit  inverse  of  the  j_Bj  matrix.  By  appropriate 
rearrangement  of  rows  and  columns  and  through  partitioning  it  can 
be  demonstrated  that  the  singularity  of  the  complete  matrix 
depends  on  the  singularity  of  a  certain  3x3  matrix.  This  3x3 
matrix  is  also  too  complicated  to  permit  its  explicit  inversion. 

Its  determinant  can  be  formulated,  however,  and  it  may  be  recalled 
that  a  criterion  for  singularity  is  whether  or  not  the  determinant 
is  zero.  The  algebraic  statement  of  the  determinant  in  question, 

D,  is 


where; 


ir-i) 

2fY{3o<-2)  +r^(.»<r-i)  +  (i-2o<)  +^2(1  +Y^)(i-r)] 

(r-i)  [  3r-2«^r(r+i)  -  +r+i)] 

+r^  {i-2c<D  +  2<^  -1 

2'r{i-r)  (i-c^r)[  r  (i-^r)+i-^3 

r  ^  (i-o<  r)^  ( r  -i)  (v-33) 
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Thus,  if  D  *  0,  the  [bJ  matrix  will  be  singular.  Equally 
Important  is  the  case  where  [b]  is  nearly  singular  (i.e.,  the 
terms  in  the  adjoint  of  the  3x3  matrix  are  very  large  in 
comparison  to  the  determinant),  since  this  will  also  produce  a 
meaningless  inverse  for  the[B3  matrix.  A  suggested  criterion  for 
this  condition  is  the  ratio  of  the  first  diagonal  element  of  the 
adjoint  to  the  determinant,  D.  Experience  has  indicated  that 
if  this  ratio  is  less  than  about  100,  then  it  is  reasonable  to 
expect  a  satisfactory  inverse  for  the  [b]  matrix.  The  first 
diagonal  element,  designated  can  be  determined  by 

*11'  (3^2^?3  -?4)  (3?3-2?2) 

+(3?2^^4«3^f4^  -  a'?2^^4«'3^«4®)  -•?4)  '  3'5’4) 

+  9^2  ^3  ^3^4^"3a^2^^’3^^3^4^"^^2  ^3^4^3  ^4^ 

+  9  ^2^  ^3^  H  -  13  ^3^  5  3®  ^4“  +  §3^ 

■  V  h  <3^  ^4^  +  1«  V  *  6  ?'2'‘^3  ^4  f  3"^ 4^ 
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D.  TRIANGULAR  PLATE  ELEMENT 

I .  Introduction 

The  formulation  of  the  triangular  plate  discrete  element 
described,  is  derived  from  and  mathematically  consistent  with,  the 
formulation  described  in  Reference  12.  The  addition  of  this 
particular  element,  as  a  companion  element  to  the  quadrilateral 
plate  described  previously,  serves  to  compliment  the  additional 
capability  available  for  the  analysis  of  shell  structures 
particularly  when  instability  analyses  are  to  be  performed. 

A  detailed  derivation  is  presented  for  the  force-displacement 
properties  of  an  orthotropic  triangular  thin  plate  element 
exhibiting  membrane  and  bending  behavior.  Included  in  these 
relationships  are  terms  for  stiffness,  stress,  thermal  stress  and 
incremental  stiffness. 
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II.  Development  of  Linear  Elastic  Membrane  Stiffness  Matrix 

A  matrix  statement  of  Castlgllano 's  first  theorem.  Part  I, 
applicable  to  the  derivation  of  discrete  element  force-displacement 
properties  may  be  written  as  follows: 

I  k]  “  [s'^j  [c]  [b]-^  (II-l) 

where  [k3  Is  the  desired  matrix  of  element  stiffness  coefficients, 

[b]  Is  a  matrix  in  which  the  rows  are  the  coefficients  of 
equations  for  the  corner  point  displacements  in  terms  of  the 
constants  of  the  assumed  displacement  pattern,  and  the  rows  of  [cj 
are  the  coefficients  of  these  same  constants  in  equations  which 
represent  the  derivatives  of  the  strain  energy  (expressed  as 
functions  of  the  constants  of  the  assumed  behavior  function)  with 
respect  to  the  respective  constants. 

For  the  triangular  plate  element,  we  assume  linear 
displacements,  l.e., 

u  =  a^  +  a^  X  +  a^y  (II-2) 

V  =  +  a^x  +  a^y  (II-3) 

where  a^,  ....  a^  are  the  constants  in  these  assumed  functions. 

Eqs.  II- 2  and  II-3  can  be  evaluated  at  corner  points  1,  2, 
and  3  (see  sketch  below)  to  yield  the  following  relationships: 
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/• 

^1 

^2 

^3 

II 

^1 

^2 

^^3 

V. 

> 

1  0 

1  X, 

1  x! 

0  o' 

0  0 

0  0 


yg 

^3 

0 

0 

0 


0 

0 

0 

1 

1 

1 


0 

0 

0 

0 

Xr 


0 

0 

0 

0 

72 

^3 


r  ^ 

ai 

^2 

^3 

"5 

^6 

mm 

L  J 

(II-4) 


whose  Inverse  is; 


[b]'^  s 


-  x^yp) 


y2^ 


-7 


3-2 

^3-2 

0 

0 

0 


0 

^3 

-X- 

0 

0 

0 


-yj 

X, 

C 

0 

0 

0 


0 

0 

0 


0 

0 

0 


(xgya-x^yg)  0 


-y 


3-2 


'3-2 


^3 

-X- 


0 

0 

0 

0 

-72 

X 


where  y^  g  “  7-^~72  ^3-2  “  ^3"^2* 

To  develop  the  Cc]  matrix  we  first  need  the  expression 
for  strain  energy  for  orthotropic  plane  stress.  This  can  be 
wirtten,  in  terms  of  strain  as; 
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(II-6) 


where  M  =  1 

and  6  ^  and  are  initial  strains.  The  strains  are 

X  J  xy 

obtained  from  the  strain  displacement  expressions  aS; 
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SO  that,  upon  expansion  and  substitution  of  Eqs.  II-7j  -8  and  -9 
Eq.  II-6  becomes 
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Setting  aside  all  terms  which  involve  initial  strains 
(these  will  be  treated  in  Section  II-B),  we  have 


The  square  matrix  on  the  right  is  the  matrix  W .  Utilizing 

[]bJ  from  Eq.  II-5  and  the  definition  of  from  Equation  II-IO 

in  Eq.  II-l  and  performing  the  indicated  operations  results  in 
the  element  stiffness  relationships  shown  on  the  next  page. 


III.  Development  of  Initial  Force  Terms 

In  the  case  of  initial  strains  (e.g.,  thermal,  strains, 
previously  accumulated  in  elastic  strains,  and  large  deflection 
strains)  it  is  necessary  to  determine  the  forces  corresponding 
to  the  initial  strains. 


With  assumed  displacement  patterns  (Eq.  II-2  and  -3)  the 
strain  energy  is  expressed  in  terms  of  the  constants  ^  a]  and  the 
initial  strains  as  shown  in  Eq.  II-6a.  The  differentiation  of 
Eq.  II-6a  with  respect  to  the  individual  constants  ^  a]  , 

Eqs.  Il-lOa  -  lOf,  results  in 

-  [«]  W  -  f#] 
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where  [  represents  the  initial  strain  terms  which  were  set 

aside  in  forming  Eq.  II -10,  and  is  given  by 
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(11-13) 


96 


The  corner  forces  on  the  element  are  obtained  by  multiply¬ 
ing  Eq.  11-13  t)y  ^  ^  a  J  =  2  -ij  ^  which  results  in 

W  “  [®'^]  [»}  -  [b'^]  f'^^] 

The  first  term  on  the  right  siae  of  Eq.  II-l4  yields  the 

corner  forces  due  to  displacement,  l.e.,  the  forces  which  would  be 

required  to  induce  the  corner  deformations,  u  and  v,  elastically. 


The  second  term  represents  the  initial  forces,  Thus 

('*]  ■  M’ m 


(11-15) 


Utilizing  from  Eq.  II-5»  Equation  11-13 

and  performing  the  product  yields  the  following  expression  for 
the  initial  forces: 
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For  thermal  strain  situations,  letting  T  represent  the 
average  temperature  change  of  the  element,  the  initial  strains 
are  defined  by 
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where  and  are  the  thermal  coefficients  of  expansion 

^  «y 

in  the  x  and  y  directions  respectively.  For  the  usual  case  of 

<^y.  =  , 

to  the  following 


=  cfC.  ,  the  initial  forces  given  by  Eq.  11-16  simplify 
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IV.  Stress  Equations 

The  equations  for  the  stresses  in  terms  of  the  corner 
point  displacements  and  initial  strains  will  he  derived  in  this 
section.  In  Section  II-A,.  it  was  noted  that  the  strains  in  the 
triangular  plate  element  are  constant  (Eqs.  11-7,  -8  and  -9). 
This  results  in  constant  stresses  which  can  be  expressed  by 
the  following: 


xy 


(II-18) 


The  first  term  on  the  right  side  of  Eq.  II-18  represents  the 
stresses  corresponding  to  the  corner  displacements  and  the  second 
term  denotes  the  initial  stresses. 

To  develop  the^S^y"]  matrix  and  initial  stresses,  the 
relationships  between  stresses  and  strains  will  be  formulated 
in  accordance  with  Hooke's  Law.  For  an  orthotropic  material 
the  strains  are  expressed  by 
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Solving  Eq.  11-19  for  the  stresses  give 
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The  first  terms  on  the  right  hand  side  of  Equations  11-20 
give  the  displacement  stresses.  By  replacing  the  strains  in 
these  terms  by  Eq.  II-7>  -8  and  -9>  the  displacement ‘stresses 
can  be  expressed  in  terms  of  the  constants  in  the  assumed  dis¬ 
placement  functions  by  the  following  _ 
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(11-21) 


Next  the  corner  point  displacements,  u  and  v,  are 
introduced  into  Eq.  11-21  by  utilizing,  from  Sect.  II-A,  the 
relationship 


f-i  ■  [•r‘ 


(II-4a) 


So  that  Eq.  11-21  may  be  written  as 


b]  b] 


(11-22) 


Where  W  represents  the  rectangular  matrix  on  the  right  hand 
side  of  Eq.  11-21  and  [bJ  is  defined  in  Eq.  II-5.  By 
explicitly  forming  the  product  of  these  two  matrices  as 
indicated,  Eq.  11-22  can  be  written  as 
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where 
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The  initial  stresses  are  now  determined  from  the  second 
terms  on  the  right  hand  side  of  Eq.  11-20.  The  consideration  of 
these  terms  gives 


^xyM 


(11-24) 


(The  minus  sign  is  not  included  in  the  terms  of  Eq.  11-24  since  it 
already  has  been  incorporated  in  the  stress  formulation  of  Eq.  II-18.) 

As  noted  in  Section  II-B,  for  thermal  conditions =o^T; 

•  • 

C  ^  =<A.J!  and  k"  0.  Then,  for  the  case  where  o<  =»<  = 
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the  thermal  stresses  are  given  by 
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V.  Development  of  Linear  Elastic  Stiffness  Matrix  (Bending) 

The  fc'"c e-displacement  properties  of  an  orthotropic 
triangular  plate  element  in  bending,  subjected  to  known  midplane 
forces,  are  derived  in  this  section.  The  element  is  pictured 
in  Figure  III-12.  As  in  previous  chapters,  the  stiffness  matrix 
is  derived  by  application  of  Castlgliano 's  Theorem.  In  the  case 
of  plates  in  flexure  subjected  to  midplane  forces,  however,  the 
stiffness  matrix  is  composed  of  two  parts,  i.e., 

[Kz]  =  [[•'fl  +  ["]]  (IV-1) 

Where 

j^KfJ  is  the  stiffness  matrix  for  flexure  alone. 

[^n  J  is  the  stiffness  matrix  associated  with  the  influence 

on  flexure  of  known  midplane  forces. 

The  present  section  is  concerned  with  the  derivation  of  the 
matrix;  the  J^n^  matrix  is  developed  in  the  next  section. 
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Figure  m -12  Triangular  Plate  Flexural  Element 
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As  shown  in  Section  III-A.3  of  Reference  12  the  [Kf] 
can  be  defined  as: 


matrix 


[kJ  =  ''  [Of]  [b] 


{IV-2) 


where,  as  In  previous  sections  [b]  Is  a‘  matrix  In  which  the 
rows  are  the  coefficients  of  equations  for  the  corner  point 
displacements  In  terms  of  the  constants  of  the  assumed  displace¬ 
ment  pattern,  and  the  rows  of  [c^J  are  the  coefficients  of 
these  same  constants  In  equations  which  represent  the  derivatives 
of  the  strain  energy  (expressed  as  functions  of  the  constants 
of  the  assumed  behavior  function)  with  respect  to  the 
respective  constants.  ' 

The  following  assumed  displacement  function  will  be  ' 
utilized  as  the  basis  for  this  derivation.  ‘ 
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where  a^, . a^  are  constants. 


(IV-3)  . 


The  angular  displacements  (slopes)  of  the  plate  are  given  by: 
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The  square  matrix  on  the  right  hand  side  of  Eq.  IV-6  is^ 
by  definition,  the  matrix. 

To  der^-lop  the  [c^]  matrix,  it  is  first  necessary  to 
express  the  flexural  strain  energy  (u^)  for  orthotropic  plates, 
in  terras  of  the  displacements:  (See  Eef.l2) 
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By  the  differentiation  of  Eq.  IV-3  and  the  substitution  of  the 
partial  derivatives  into  Eq.  IV-7,  the  following  expression  for 
the  strain  energy  in  terms  of  the  constants  of  the  displacement 
function  is  obtained: 

“f  °  3  r.  [  *  24  a^agX  +  36  ag^x®)  +  Dy  (4832 

+  5Sa„^y^  +  4a„^x^  +  24aoa„y  +  24a„aQ  +  8a,aQx) 

f  (  if  r  y  3  y  (tv.8) 

+  S^yx^^x  +  12a^a.^y  +  4a^a^x  +  12a2agX 

2  *1 

+  36aga^xy  +  12aga^x  )  +  2D^y  (ag^  +  4aga^y  +  4a^^y^)J  dA 
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VI,  Development  of  Incremental  Stiffness  Matrix 


The  [^nj  matrix  can  also  be  derived  through  application  of 
Castigliano 's  First  Theorem,  Part  I  (See  Reference  12) .  The  procedure 
is  represented  by  the  relationship 


[n]  =  [b-1]  ''  [o„]  [b] 


The  strain  energy  expression  to  be  employed  in  the  formulation 
of  the  matrix  is; 
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(IV-12) 

Where  N  ,  N  and  N  are  the  known  midplane  forces,  these  forces 

^  Jr 

will  have  been  evaluated  by  performance  of  an  independent  midplane 
displacement  analysis,  wherein  the  element  is  assumed  to  sustain 
constant  midplane  stresses  (  =  a^,  cr^  =  =  a^) .  Since 

N  =  her  N  =  her  and  N  =  h.T  ,  the  midplane  forces  are  also 

A  A  y  jf  Ay  Ay 

constant  by  taking  note  of  this  consideration,  utilizing  the  displacement 

function  of  Eq.  IV-2,  and  performing  the  operations  indicated  by 

Eq.  IV-12,  one  obtained;  U„  =  U  +  U  +  U 
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Considering  each  of  the  energy  components  separately,  the 
partial  derivatives  of  the  energy  with  respect  to  the  constants 
may  be  stated  as  follows; 
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The  three  rc„l  matrices  are  shown  on  the  following  pages.  The 
i  1  ^ 

T  w  vi*t  v^rr  4  v^  *n  a  4*  v»4  4  4>l>k«xk 


X  y 


terms  appearing  in  these  matrices  are  defined  in  the 


following  section  -  Geometric  Properties. 

The  I^C^j  matrix  in  Eq.  IV-16  equals  the  sum  of  the 
three [c^]  matrices  of  Eqs.  XV-l4a,  l4b  and  l4c.  Hence,  the 
incremental  stiffness  matrix  [n]  is  experssed  as; 


-1 


(IV-lla) 
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VII.  Geometric  Properties 

The  [c^J  and  matrices  contain  a  group  of  terras 

resulting  from  integration  of  the  energy  expressions.  These 
terms  are  defined  as; 

"  Ja  ^  (IV-15) 

Thus,  they  are  simply  geometric  properties  of  the  element.  Many 
of  them  are  well  known  section  properties. 

In  explicitly  formulating  the  terms  it  was  convenient 

to  use  previously  determined  ^  terms  with  reference  to  the 

coordinates,  as  shown  in  the  following  sketch,  and  transform 
these  to  the  x-y  system. 


The  coordinate  transformation  from  the  x-y  system  into  the 
coordinates  for  any  given  point  is 


cos  0 
-sin  0 


sin  0 
cos  0 


(IV-16) 


where  sin  0  = 


cos  0 


Alternatively  the  coordinate  transformation  from  the  system 

to  the  x-y  system  is  given  by 


(IV-l6a) 
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All  the  remaining  I  „^'s  are  determined  by  a  similar  procedure. 

X  y  <  ,• 

A  listing  of  the  necessary  y'^'s  is  given  on  Fig.iii-13. 

The  area.  A,  and  moment  properties  about  the  ^  and  ^  axes 
were  determined  by  direct  integration  within  the  proper  limits 
to  yield  the  expressions  shown  on  Fig.  III-14. 
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VIII,  Development  of  Corner  Thermal  Moments. 

Reference  l4  has  established  a  procedure  for  the  derivation 
of  thermal  forces  that  is  consistent  with  the  procedures  employed 
in  deriving  the  and  £nj  matrices,  i.e.,  a  procedure  based 
on  Castigliano's  First  Theorem,  it  is  simpler  and  appears  equally 
accurate,  however,  to  derive  the  thermal  forces  by  means  of .  the 
scheme  to  be  described  in  the  following.  •  '  ' 

In  developing  the  corner  thermal  moments,  the  average  distri¬ 
buted  thermal  moments  (about  the  x  axis)  and'  My*=^  .  (about 

the  y  6ixis)  are  taken  as  constant  throughput  the  element.  The 
"lumped"  corner  thermal  moments  are  based  upon  these  average 
moments.  The  average  moments  are  defined  as  the  arithmetic 
average  of  the  distributed  moments  at  the.  three  corners ^ of  the 
triangular  element.  Thus 


+  Myg-" 


,(IV-18) 


where  and  are  the  distributed  moments  at  the  corners 


due  to  the  temperature  gradient,  through  the  thickness  of. the 
plate.  These  distributed  corner  moments  are  defined  as: 
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^  is  a  thickness  coordinate  measured  positively  in  the  positive  z 
direction  from  the  neutral  axis  of  the  plate. 
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It  is  assuhied  that  the  average  thermal  moments  are  distributed 
around  the  edges  of  an  arbitrary  triangle  as  shown  below; 


«  M 

^  X 


The  , distributed  thermal  moments  are  concentrated  at  the  corners 
of  the  triangle  byassignir  ■  one  half  of  the  total  moment  along 
a  side  tp  each  of  the  corners  bounding  the  side.  For  example: 
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The  thermal  corner  forces  are  zero. so  that  the  thermal  moments 

2  .  I 

and  forces  are  expressed  in  matrix  notation  by 
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IX.  Development  of  Stress  Resultant  Bending  and  Twisting 
Moments  and  Vertical  Shear  Forces 

It  appears  most  efficient,  from  the  standpoint  of  practical 

interpretation,  to  express  the  "stresses"  as  moments  and  shears 

per  unit  length.  Hence,  the  bending  moments,  M  and  M  ,  the 

X  y 

twisting  moment  and  the  shear  forces,  and  Q^,  shown  iu 

the  sketch  can  each  be  calculated  during  an  analysis.  They  are  to 
be  computed  at  the  centroid  of  the  element.  The  stresses,  which 
are  dependent  upon  the  construction  details  of  the  cross-section, 
can  then  be  hand  calculated  from  the  moments  and  shear  forces. 

The  moments  and  forces  may  be  expressed  in  terms  of  the 
deflected  surface,  as; 
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where  D  ,  D  and  D  are  defined  below  Eq.  IV-7  and 
X  y  xy 

D-  =  D  >«'  +  D 

Q  X  ^yx  xy 

By  differentiating  the  assumed  displacement  function  (Eq.  IV-3), 
the  following  derivatives  are  obtained 
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(IV-23) 


ag  +  2  a^y 
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Substitucing  Eqs.  IV-23  into  Eqs.  IV-22  and  evaluating 
the  expressions  at  the  centroidal  coordinates  (x,  y),  yields  the 
following  matrix  equation  for  the  moments  and  forces: 
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(IV-24) 


where: 

[Sz]=  - 


and 


0 

2D 

y 

0  2y**  D 

^y  y 

^'liyV 

6Dyy 

0 

2DyX 

0 

2/^  D 
yx  X 

0 

2Dx 

en^s 

6y^  D  V 
'^x  x*^ 

0 

0 

0 

0 

0 

0 

0 

"^xy 

■^^xy^ 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

6d 

y 

0 

0 

X  = 

+  x^) 

y  = 

3  (j^2 

+  y3) 

From  Eq.  IV-6  it  is  noted  that; 


(IV-6a) 
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SO  that  Eq.  IV-24  may  he  written  as 
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xy 
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(IV- 24a) 


which  permits  the  moments  and  forces  to  he  calculated  from  the 
previously  computecf  j;  is  placements.  To  account  for  the  presence  of 
thermal  moments,  Eq.  IV-24a  is  modified  as  follows: 
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Where  [s,]  =  [g^]  [b]'^  and 
thermal  moments  and  forces,  is  given  by: 


the  column  of 


#  I 

In  Eq.  IV-26,  M  and  are  the  average  distributed  thermal 

y 

moments  defined  by  Equation  IV-I8. 
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f.xV  =  +  Ijl*)  ®^“^®  =“^6  +  2  (I^3j  -  1.^3) 

(sin  e  cos^e  -  Sln^e  cos  0)  +  I  „  .  (cos^0  -  4  sin®0  cos^0 
+  Sin  0) 


^x  y  ^^;5»4  “  31;^  2  ^2'  ®  cos^0  +  3  (i^^3  -  ^^35^ 

sin^0  cos^Q  +  (3  2^2  “  4)  ® 

4  4 

-  o  sin  0+1  ,  cos  0 

^^4' 


3  _ 


xy 


(I;5,4  ‘  31^2^2)  ®  ® 


>? 

+  3  (I 


7  3^  - 


+  (31^  2^2  "  ^$4)  ®  cos^0 

li  ii 

-  I  .  Q  sin  0  +  I  .  _  cos  0 

f§3 


FIGURE  m- 13  (CONTINUED) 
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E.  FRAME  ELEMENT 

I.  Introduction 

The  formulation  of  the  additional  frame  discrete  element 
which  has  been  incorporated  into  the  MAGIC  II  System  is 
essentially  Identical  to  that  described  in  the  original  MAGIC 
Engineer's  Manual  (Reference  1).  All  element  matrices  available 
to  the  original  frame  element  are  available  to  this  frame  as 
well,  i.e..  Stiffness,  Distributed  Pressure,  Thermal  Load, 
and  Consistent  Mass. 

The  addition  of  this  element  is  primarily  intended  to 
serve  the  purpose  of  providing  a  companion  frame  element  to  the 
quadrilateral  and  triangular  plate  elements  which  have  been 
added  to  MAGIC  II.  The  use  of  this  element  in  conjunction  with 
the  newly  added  quadrilateral  and  triangular  plate  elements 
will  provide  a  powerful  capability  for  linear  eigenvalue 
stability  analyses  of  stiffened  shell  structures. 

II.  Additional  Element  Matrices 

As  pointed  out  previously,  all  element  matrices  available 
to  the  frame  are  described  in  detail  in  Reference  1.  An 
additional  incremental  stiffness  matrix  has  been  provided  and  is 
presented  here  (References  l4  and  15) . 

For  the  frame  element  (Fig.  Ill -15)  linear  polynomial  axial 
and  torsional  displacement  mode  shapes  are  constructed  while  a 
cubic  polynomial  displacement  mode  shape  is  constructed  for 
flexure  in  each  of  the  two  principal  planes  of  bending.  The 
above  mentioned  mode  shapes  are  assumed  to  take  the  following  form 


(2.1) 

(2.2) 

(2.3) 

(2.4) 


The  above  mode  shapes  lead  to  a  total  of  12  undetermined 
coefficients  for  the  element  which  are  chosen  to  correspond 
to  three  translational  and  three  rotational  displacement  degrees 
of  freedom  at  each  end  of  the  element.  A  transformation  from 
generalized  coordinates  to  grid  point  displacement  degrees  of 
freedom  is  effected  by  writing  (at  X  =  0) 
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and  at  X  =  L 
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It  is  to  be  noted  that  the  are  referred  to  as  field  coor¬ 

dinate  displacement  degroes  of  freedom. 

Upon  analytical  inversion  of  equation  (2.7)  we  have  the 
desired  relationship  between  the  and  jcTj  displacement  vectors. 

“  r/?yJ  {<^]  (2.10) 
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The  strain-displacement  relationships  for  the  frame 
element  can  be  written  in  the  following  form: 

£  .  (u^  +  1/2  v/  -  +  1/2  w^2  _  2W^^)  (2.11) 


The  total  potential  energy  functional,  J  ,  which  arises 
in  consequence  of  the  strain  relations  of  equation  (2.11)  is 
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The  first  integral  in  Equation  (2.12)  consists  of  the  well 
known  linear  membrane  and  flexure  terms  respectively  while  the 
second  integral  arises  from  the  retention  of  the  quadratic  terms 
in  the  strain-displacement  relation. 

In  consideration  of  the  non-linear  portion  of  Equation  (2.12) 
it  is  noted  that  the  first  term  is  the  non-linear  membrane-flexure 
coupling  term.  This  contribution  gives  rise  to  terms  which 
adversely  affect  the  element  linear  membrane  and  flexure  stiffness 
and  is  the  term  which  will  be  considered. 

If  consideration  is  given  to  work  done  by  midplane 
loads  during  the  displacement  of  the  structure  during  lateral 
loads,  the  following  may  be  written  for  the  first  non-li,iear 
contribution  of  Equation  (2.12) 


I'  = 


/  (''x'  -  “x") 


dx 


(2.13) 


where  P  is  the  axial  stress  resultant  which  is  a  known  quantity. 
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From  the  assumed  displacement  functions  (Equations  2.2 
and  2.3)  differentiate  and  obtain: 


=  b^  +  2b2X  +  Sb^x 

r 

+  2C2X  +  Sc^x' 


(2.14) 


substitution  into  Equation  (2.13)  and  integrating  obtain: 
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+  3b2l3  3L^  +  9/5  b^^L^  +  c^^L  +  2c^C2L' 
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Upon  taking  the  partial  derivatives  of  J  with 
respect  to  the  coefficients  b^  thru  b^  and  c^  thru  Cy  the 

incremental  stiffness  matrix  [A/]  is  obtained  referenced  to 
field  coordinate  displacement  degrees  of  freedom  and  is  shown  in 
Equation  (2.16). 


Noting  from  Equation  (2.10)  that 


r-v2j 


is  obtained  as  follows: 

L^^iJ -- [/Z^] 
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(2.17) 


when  \_/V  1 ]  is  referenced  to  grid  point  displacement 

degrees  of  freedom. 


SECTION  IV 

COMPUTATIONAL  PROCEDURES 


A.  Introduction 

The  MAGIC  II  System  for  Structural  Analysis  offers  a 
variety  of  computational  procedures  to  the  User.  Among  these 
are  the  capability  to  perform  static  analyses,  statics  with 
condensation,  statics  with  prescribed  displacements,  stability, 
dj-namics  (modes  and  frequencies)  and  dynamics  with  condensation. 

The  proper  usage  of  these  procedures  in  the  context  of  performing 
actual  structural  analyses  is  described  in  detail  in  Volume  II 
of  this  document  (The  User's  Manual). 

In  addition,  the  powerful  matrix  abstraction  capability 
built  into  the  MAGIC  II  System  allows  analyses  to  be  performed  which 
require  the  use  of  Static  and  Dynamic  Substructuring.  In  order 
to  clarify  the  operations  involved  utilizing  these  approaches, 
a  detailed  presentation  follows. 

B.  Static  Substructuring 

A  primary  attraction  of  the  matrix  methods  of  structural 
analysis  is  that  many  significant  problems  can  be  solved  with  very 
modest  computer  programs.  In  the  simplest  case,  a  capability  to 
generate  and  assemble  a  certain  type  of  element  stiffness  matrix  and 
solve  the  resulting  system  stiffness  equation  is  sufficient.  On 
the  other  hand,  the  automated  analysis  systems  required  to  cope  with 
the  large  classes  of  structures  in  a  practical  design  situation  bear 
little  relation  to  the  type  of  finite  element  computer  program 
mentioned  above.  Practical  analysis  tools  must  be  Impleraented  as 
an  integral  part  of  the  overall  structural  analysis  and  design 
cycle. 
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It  is  apparent  from  the  size  and  complexity  of  many 
structures  that  a  great  deal  of  data  is  involved-  The  detailed 
design  specification  of  the  structure  is  spread  over  numerous 
drawings  and  through  many  documents.  A  realistic  physical  model 
necessitates  thousands  of  gridpoints  and  finite  elements. 

Translated  into  computer  program  input,  related  data  items 
include  coordinates  of  all  gridpoints,  degree-of-freedom  specifica¬ 
tions  of  all  gridpoints,  connection  specifications  for  all  finite 
elements,  etc.  This  volume  of  input  data  implies  a  need  for  many 
hours  of  computer  time  for  numerical  solution.  Finally,  extensive 
output  inevitably  results  from  such  an  analysis. 

Effective  management  of  the  voluminous  data  associated 
with  structures  of  this  type  is  usually  the  decisive  consideration  in 
establishing  the  basic  analysis  process.  For  these  reasons  an 
analysis  can  be  undertaken  by  substructuring  (References  16,  17). 

In  general,  the  substructuring  process  proceeds  in  four  major 
phases  as  outlined  below. 

Phase  I  is  concerned  with  the  individual  substructures  of 
the  total  Structure.  In  the  Phase  I  analysis,  each  substructure 
is  considered  individually.  Input  data  is  prepared  and  calculation 
is  carried  forward  to  determine  matrix  representations  referenced 
to  the  substructure  interfaces. 

Phase  II  considers  the  structure  as  a  whole.  The  interface 
stiffness  matrices  for  the  individual  substructures  are  assembled 
and  the  complete  set  of  interface  displacements  is  determined. 

Based  upon  interface  displacements  obtained  from  Phase  II 
and  auxiliary  Information  from  the  individual  Phase  I  analysis. 

Phase  III  completes  the  conventional  finite  element  analysis.  Each 
substructure  is  considered  in  turn.  Prediction  of  the  primary 
displacement  variables  is  completed  and  secondary  variables  such  as 
forces  and  stresses  are  calculated. 
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Phase  IV  is  a  nonintegral  step  designed  to  translate  the 
conventional  Phase  III  results  into  a  form  desired  hy  the  stress 
analyst  for  the  determination  of  margins  of  safety. 

1.  Phase  I 

The  matrix  algebra  of  the  Phase  I  analysis  is  deceptively 
straight  forward  when  reduced  to  the  essential  calculations 
relevant  to  the  primary  displacement  variables  and  stripped  of  the 
problematical  data  storage  and  retrieval  steps  inherent  in  the 
computer  program.  The  simplified  symbolic  statement  of  the  process 
is  considered  appropriate  in  the  present  context.  Accepting  this 
viewpoint,  the  first  step  of  the  Phase  I  analysis  process  yields 
potential  energy  expressions  for  the  individual  finite  elements 
given  by, 

V  ■  jLs.jck,]  {6,}  -  (p^) 

where 

[Kg]  is  the  element  stiffness  matrix, 

{6g}  is  the  relevant  gridpoint  displacement  vector,  and 
{P^}  is  the  element  total  applied  load  vector. 

These  individual  element  potential  energies  are  then  assembled 
to  form  a  potential  energy  expression  for  the  substructure  under 
consideration. 

♦p  •  I  L6J  [K]{6}  -  L6J  {P>  (2-2) 

where 

[K]  is  the  substructure  stiffness  matrix, 

{6}  is  the  substructure  displacement  vector,  and 

{p}  is  the  total  substructure  applied  load. 
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The  Phase  I  analysis  is  carried  forward  by  rewriting  the  sub¬ 
structure  potential  energy  in  partitioned  form  to  refl-ect  the 
division  between  interior  gridpoint  degrees -of -freedom  W 
and  interface  (boundary)  gridpoint  degrees -of -f reedom  [</^j  i.e.. 


Contributions  to  the  potential  energy  vAiich  stem  from  the 
interior  gridpoints  are  complete  while  additional  contributions 
will  be  added  in  at  the  interface  gridpoints  upon  assembly  of  the 
substructures.  Adv.'ntage  is  taken  of  the  completeness  at  the  interior 
points  by  solving  for  these  displacements  degrees -of -freedom  in 
terms  of  the  interface  degrees -of -freedom.  The  result  is  given  by. 


{6^}  » 


(2-4) 


Backsubstitution  of  this  relation  into  Equation  2-3  yields  the 
objective  substructure  potential  energy  expression  referenced  in 
interface  degrees -of -freedom,  i.e. 


where 
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[Kjj]  is  the  substructure  interface  stiffness  matrix, 

{6  }  is  the  substructure  interface  displacement  vector,  and 

D 

{Pw}  is  the  substructure  interface  load  vector, 
b 
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The  individual  substructure  potential  energy  expressions  of  the 
form  defined  in  Eauation  2-5  are  the  basic  Phase  I  analysis  results 
required  to  construct  the  governing  stiffness  equation  for  the 
entire  structure  in  the  Phase  II  analysis  process. 

The  foregoing  statement  of  the  Phase  I  analysis  process 
actually  implies  a  complete  general  purpose  computer  program  for 
stress  analysis  plus  the  capability  to  form  and  store  such  items 
as  interface  stiffness  matrices  on  magnetic  tape  for  subsequent 
access .  It  is  instructive  to  take  the  viewpoint  of  the  structural 
analyst  and  re-examine  the  Phase  I  analysis  process  as  an  applica¬ 
tion  of  the  MAGIC  II  Structural  Analysis  System. 

By  definition.  Phase  I  proceeds  against  a  number  of  sub¬ 
structures  of  the  total  structure  although  inclusion  of  the 
complete  structure  within  a  single  substructure  would  yield  a 
conventional  one-pass  linear  stress  analysis.  The  reasons  for 
division  of  a  structure  into  multiple  substructures  are  many  and 
varied.  Unnecessary  breakdown  into  substructures  is,  of  course, 
inefficient. 

A  primary  reason  for  substructuring  stems  from  the  fact  that  it 
is  efficient  to  confirm  a  large  quantity  of  input  data  via  subsets. 

With  substructuring,  an  analyst  can  forcus  his  attention  upon  a 
limited  region  or  component  in  specifying  input  data.  This  subset 
of  data  can  then  be  processed  through  data  checking  executions. 

Such  executions  involve  only  the  relatively  small  quantities  of  data 
of  current  interest  with  the  result  that  turn-around  is  rapid  and 
inexpensive. 

Substructuring  can  also  shorten  the  calendar  time  required 
to  confirm  the  input  data  for  a  large  structure.  The  roason  for 
this  is  that  substructuring  facilitates  distribution  of  the  input 
data  specification  effort  to  a  number  of  analysts  for  nearly  independent 
simultaneous  preparation.  It  is  worth  mentioning  that  the  automatic 
generation  of  structural  plots  fro;u  the  input  data  is  an  important 
aid  to  the  confirmation  of  input.  Plots  of  the  structural  components 
taken  individually  are  desirable. 
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The  benefits  of  substructuring  large  scale  structures 
extend  beyond  the  input  data  confirmation  stage  through  execution. 
The  effective  matrix  banding  may  be  improved  by  substructuring. 

Long  continuous  executions  are  avoided.  Numerous  restart  points 
are  automatically  provided.  Most  important,  the  Phase  I  executions 
are  spread  over  the  period  of  time  required  to  complete  the 
specification  of  data  for  all  the  component  substructures. 
Executions  in  the  succeeding  phases  may  be  similarly  spread  to 
generate  results  paced  by  progress  in  evaluation. 

Subs true turing  is  particularly  advantageous  when  localized 
modifications  in  structure  or  applied  loading  arise  subsequent  to 
the  analysis.  Such  modifications  can  often  be  accommodated  by 
re-analysis  of  only  those  substructures  affected. 

Having  discussed  the  motivation  for  and  the  scope  of  the 
Phase  I  analysis,  the  following  paragraphs  focus  upon  the  several 
steps  involved.  Preprinted  input  forms  are  employed  to  simplify 
the  specification  of  input  data.  These  forms  are  designed  to 
provide  automatic  internal  generation  of  data  whenever  possible. 

For  example,  repetitious  data  need  only  to  be  specified  initially 
followed  by  any  exceptions.  (See  Volume  II  -  User’s  Manual). 

The  first  executions  of  the  MAGIC  II  Analysis  System  are 
undertaken  to  confirm  the  input  data  deck  as  discussed  earlier. 

The  deck  is  read  and  the  Implied  data  is  generated  explicitly. 
Consistency  of  the  data  is  checked  and  all  data  items  are  stored 
for  execution  restart  and  printed  for  further  checking  by  the 
analyst.  In  addition,  a  magnetic  tape  is  generated  for  automatic 
plotting  of  the  finite  element  model. 

Upon  acceptance  of  the  input  data  specification  by  the 
analyst,  the  actual  Phase  I  analysis  is  undertaken  for  the 
substructure  under  consideration.  This  analysis  is  a  complete 
linear  stress  analysis  of  the  substructure  under  the  assumption 
that  the  interface  displacements  are  completely  fixed.  The  output 
obtained  from  this  analysis  provides  further  important  confirmation 
of  the  finite  element  model.  Moreover,  these  results  often  provide 
useful  preliminary  information  about  substructure  behavior. 


In  addition  to  the  preliminary  stress  analysis  results,  the 
Phase  I  analysis  generates  and  stores  the  interface  referenced 
stiffness  and  applied  load  matrices  as  well  as  the  other 
information  required  in  subsequent  analysis  phases. 


2.  Phase  II 


The  Phase  II  analysis  begins  with  the  substructure  Interface 
matrices  from  Phase  I  and  carries  the  analysis  process  through 
prediction  of  the  interface  displacement  variables.  Phase  II  is 
the  only  part  of  the  analysis  process  which  deals  with  more  than 
one  substructure  at  a  time. 


The  substructure  potential  energy  expressions  (Equation 
2-5)  are  the  point  of  departure.  Such  an  expression  is  known 
for  each  substructure,  e.g. 


=  1 
pb  2 


n 


(2-8) 


The  interface  displacements  pertinent  to 

are  known  as  elements  of  the  total  list 
{a]  for  the  assembled  structure.  This 
expressible  mathematically  by  a  Boolean 

each  substructure 

of  interface  displacements 
relationship  is 

transformation . 

Symbolically, 

(2-9) 

Introduction  of  this  transformation  into  the  set  of  independent 
interface  displacement  degrees -of -freedom  for  the  total  structure  yields 

*pb'  = 

j=l»2,  .  .  .  n  (2-10) 

where 

cl  D  & 

(2-11) 

(2-12) 
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The  substructure  interface  potential  energy  expressions  of  the 
form  of  Equation  2-10  are  added  to  obtain  the  complete  potential 
energy  for  the  structure,  i.e.. 


where 


$ 

p 

'  [k]{a}  -  {p} 

(2-13) 

CK]  = 

(2-14) 

J 

{P} 

2{p(J )  J 

(2-15) 

J 


The  formality  of  transforming  to  conformable  substructure 
matrices  before  assembly  of  the  total  structure  matrices  is 
avoided  in  actual  practice.  Instead,  a  nonconformable  sum  is 
effected  to  obtain  Equation  2-13  directly  from  Equation  2-8. 

The  objective  equation  governing  displacement  of  the 
interface  points  follows  Immediately  by  executing  the  variation 
rf  the  potential  energy  expression  cf  Equation  2-13.  The 
result,  retaining  the  symbolism  of  a  single  load  condition, 
is  given  by, 

[k]  Ca)  =  (p)  (2-16) 

The  total  interface  matrix  [k3  is  a  symmetric  matrix  which 
is  stored  in  banded  form.  Solution  is  effected  by  triangulari- 
zation  of  the  matrix  and  back  substitution  for  the  set  of 
interface  displacement  variables  appropriate  to  each  load 
condition. 


3.  Phase  III 

Phase  III  picks  up  the  matrix  descriptions  of  the  indivi¬ 
dual  substructures  generated  in  Phase  I  and  the  solution  for  the 
interface  displacements  from  Phase  II  and  carries  the  analysis 
process  forward.  Firstly,  the  relevant  interface  displace¬ 
ments  are  extracted  from  the  complete  set  (Equation  2-9).  Then, 
the  interior  displacements  are  calculated  using  Equation  2-4, 
i.e. , 

{6i}  =  (2-17) 

With  this  result  all  the  primary  variables  for  a  given  sub¬ 
structure  are  known  and  various  secondary  items  are  computed. 

For  example,  stresses  are  available  for  each  finite  element  in 
the  substructure  via  a  relation  of  the  form 

{o}  =  [S]{6}  -  {^}  (2-18) 

where 

(o)  is  the  element  stress  vector, 

[S]  is  the  element  stress  matrix,  and 

i4>}  is  the  thermal  stress  correction  vector. 

Many  useful  additional  items  are  calculated  in  Phase  III 
in  the  MAGIC  II  System,  Included  are:  element  forces,  reactions 
and  force  balance.  Even  so,  this  information  falls  short  of  that 
desired  for  the  margin  of  safety  determinations.  This  gap 
between  the  conventional  finite  element  stress  analysis  results 
and  the  information  desired  by  the  stress  analyst  necessitated 
the  extension  of  the  automated  analysis  process  to  include  a 
Phase  IV. 


I4l 


4. 


Phase  IV 


Phase  III  was  set  up  to  compute  the  normal  finite 
element  results  for  substructures  specified  by  the  stress 
analyst.  These  results  are  stored  on  magnetic  tape  as  well 
as  printed.  This  magnetic  tape  furnishes  the  primary  input 
data  for  the  Phase  IV  analysis.  Phase  IV  is  initiated  either 
following  Phase  III  or  after  examination  of  the  printed  results 
from  Phase  III  at  the  discretion  of  the  stress  analyst. 

The  computations  of  Phase  IV  are  designed  to  automatically 
reduce  the  predicted  behavior  data  for  evaluation  by  the  stress 
analyst  in  terms  of  margins  of  safety.  A  typical  computation  of 
this  automatic  reduction  process  is  the  consideration  of  tensile 
yielding.  This  involves  the  interpolation  of  the  allowable 
stress  from  the  appropriate  temperature  referenced  table  on 
a  magnetic  tape  file.  Then,  the  equivalent  stress  is  calculated 
from  the  actual  multi-axial  stress  state.  This  equivalent 
stress  state  is  interpreted  via  Von  Mises  yield  criterion  for 
comparison  with  the  allowable  stress.  Comparison  is  made 
quantitatively  in  terms  of  a  margin  of  safety.  The  results 
of  this  comparison  are  printed  with  explicit  labelling  and 
an  asterisk  (*)  is  employed  to  identify  all  negative  margins 
of  safety.  At  the  present  time.  Phase  IV  is  a  non-integral 
part  of  the  MAGIC  II  System. 
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C.  Dynamic  Substructuf Ing 

The  development  of  the  concept  of  dynamic  suhstructuring^^®^ 
which  includes  means  for  reducing  the  number  of  degrees  of  freedom 
in  a  structural  dynamics  system,  is  presented  in  this  section. 

Nine  distinct  phases  ranging  from  discussion  of  finite  element 
building  blocks  in  Phase  1  to  computation  of  system  modes  and 
frequencies  in  phase  9  are  defined. 

The  powerful  matrix  abstraction  capability  built  into 
the  MAGIC  II  System  makes  possible  the  employment  of  the 
computational  procedure  which  is  outlined  below. 

PHASE  1  -  ENERGY  FUNCTIONALS 

The  basic  building  blocks  of  the  mathematical  model  for  a 
complete  structural  system  are  taken  to  be  the  finite  element 
energy  functions.  The  potential  energy  functional  for  the 
continuum  of  an  individual  finite  element  is  discretized  by  the 
construction  of  assumed  modes  in  accord  with  the  Ritz  procedure 
Presuming  an  admissible  assumed  mode  discretization  of  the  k^" 
finite  element  of  the  form 

=  (1) 

The  resulting  algebraic  expressions  for  the  element  energy  functions 
based  on  the  concept  of  consistent  matrices  for  stiffness  ’ 

applied  loading  ^ F^ j  ,  damping  and  mass  may  be  cast 

into  matrix  notation  as  follows : 

a.  Strain  Energy 
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b.  External  Work 


(3) 


c.  Dissipation  Energy 


or 


^^(k)  ^1.  |^^^(k)j  j-  ^(wjj  8  jk)J 

°T  (stractural  damping) 


(4) 


(5) 


d.  Kinetic  Energy 


4  (k)  1  i8,k)ir„(k)''/g  (k)| 

M  2LeJLeJLeJ 


(6) 


The  equations  governing  dynamic  response  of  a  structural  system  are  derivable 

(201 

from  such  energy  functions  via  a  generalized  Lagrange  equation  .  Total  energy  func- 
tious  for  the  complete  structural  system  are  required  in  this  formulative  process.  Given 
these,. the  objective  matrix  equation  is  readily  obtained  as 


(7) 


Nine  analysis  phases  trace  the  dynamic  substructuring  process  from  the  basic  finite  ele¬ 
ment  energy  functions,  through  the  construction  of  total  energy  functions  for  the  complete 
structural  system,  to  a  governing  matrix  equation  of  the  form  of  Equation  7,  which  is 
particularly  well  suited  to  dynamic  response  analyses. 

Development  toward  the  objective  expressions  of  the  several  total  energy 
functions  for  the  complete  structural  system  is  carried  forward  in  each  of  the  subsequent 
analysis  phases  by  the  construction  of  a  displacement  coordinate  transformation  of  the 
form 


(8) 
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This  transformation  relation  is  employed  to  change  the  energy  functions  from  expression 
in  terms  of  the  generalized  coordinates  {h)  to  expression  in  terms  of  a  new  set  of 
generalized  coord Imites 

It  is  useful  to  illustrate  the  impact  of  such  a  transformation  on  a  representative 
energy  form  before  proceeding.  Beginning  from  an  energy  expression  given  by 

Introduction  of  a  transformation  of  the  form  of  Equation  8  yields 


The  first  term  (a)  in  this  result  can  be  associated  with  the  reference  level  of 
the  energy  and  discarded.  The  second  term  (b)  defines  the  modified  nucleus  matrix 
occasioned  by  the  coord^<nate  transformation.  Finally,  the  third  term  (c)  defines  a  gen¬ 
eralized  work  contribution  that  arises  out  of  the  coordinate  transformation.  Only  the 
second  term  is  written  explicitly  in  denoting  such  transformations  in  the  text  of  subsequent 
analysis  phases,  i.e. 


where 


The  change  in  reference  levf '  and  the  contribution  to  the  work  energy  are  assumed.  The 
above  explanation  permits  conciseness  in  subsequent  steps  without  ambiguity. 
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PHASE  2  -  FINITE  ELEMENT  ASSEMBLY  SUBSTRUCTURE  LEVEL 


This  portion  of  the  analysis  concerns  the  assembly  of  the  individual  finite  element 
energy  functions  that  comprise  a  typical  substructure.  The  element  gridpoint  displacement 
degrees  of  freedom  |  ^  *^re  assumed  to  be  referenced  to  compatible  reference  axes. 
Under  this  assumption  the  displacement  set  for  the  substructure  |  B  j.  is  related 
to  that  of  its  element  a  Boolean  transformation  of  the  form 


{s  [r  «]{s  (i)}. 


(13) 


th 


Introduction  of  this  coordinate  transformation  into  Equation  2  yields  the  k  element 


strain  energy  with  reference  to  the  complete  set  of  gridpoint  displacement  degrees  of 


.th 


freedom  for  the  j  substructure,  i.e. 


(14) 


Summation  over  all  of  the  finite  elements  yields  the  strain  energy  expression  for 
the  assembled  substructure  as 


(15) 


where 


(16) 


All  of  the  substructure  energy  forms  are  derived  similarly  by  introduction  of  the  coordinate 
transformation  of  Equation  8.  Explicit  statement  of  the  matrix  algebra  for  each  of  the 
energy  functions  is  omitted  to  avoid  needless  repetition. 

The  consistent  structural  mass  matrix  of  the  finite  element  model  usually  requires 
modification  to  incorporate  nonstructural  masses  attached  to  the  structure.  Accordingly, 
provision  is  made  in  this  Phase  2  analysis  to  check  and  augment  the  substructure  consistent 
structural  mass  matrix. 
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Displacements  that  are  prescribed  as  functions  of  time  can  be  eliminated  as  degprees 
of  freedom  at  this  point  and  carried  forward  as  generalized  applied  loads.  The  preparation 
for  this  elimination  is  to  rearrangi>  the  set  of  gridpoint  displacements  to  exclude  displace¬ 
ments  that  are  prescribed  zero  (fixed  displacements)  by  setting  corresponding  matrix  rows 
and  columns  to  zero  and  position  those  prescribed  nonzero  below  the  retained 

degrees  oi  freedom  |  S  .  The  conformably  partitioned  strain  energy  is  given  by 

♦u'2  [l*iaj-l*ibjj  [K,]  jh.b]  f{8la} 

- I - J - -  (17) 

{V} 


The  reduction  occasioned  by  prescribed  displacements  is  approached  via  the  general 
transformation  of  Equation  8.  In  the  interest  of  uniformity,  we  write 

where 

{82}'{8ia}  <“> 

Substitution  of  this  prescribed  displacement  transformation  into  the  energy  form  of 
Equation  17  yields  a  modified  quadratic  form  given  by 

♦u4i*2J[s]{®2} 


where 
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The  associated  contribution  to  the  work  arises  as 

*W=  U2J{>’2} 

where 

{>’2}  =  [r2]'"[''i]{’'2}  ■  h,b]{Sib}- 

It  should  be  noted  that  generality  may  be  lost  unnecessarily  by  invoking  the  pre¬ 
scribed  displacement  reduction  at  this  point.  Generally,  only  relatively  small  reduction 
in  the  order  of  the  matrices  of  the  probl'im  is  realized.  It  is  best  to  defer  such  speciali¬ 
zations  of  the  model  as  long  as  possible  to  preserve  its  generality.  The  procedure  outlined 
here  to  effect  the  prescribed  displacement  reduction  is  applicable  at  whatever  stage  the 
reduction  is  carried  out. 

PHASE  3  -  CONDENSATION  (SUBSTRUCTURE  LEVEL) 

This  phase  of  the  analysis  process  derives  from  the  likelihood  that  the  complete 
set  of  gridpoint  displacement  degrees  of  freedom  -  {82}  are  not  essential  to  the  objective 
structural  dynamics  analyses.  For  example,  the  gridpoints  in  the  finite  element  model 
may  have  been  dictated  by  the  natural  breakdown  of  the  structure  into  components,  or  the 
intended  use  of  the  model  for  stress  analyses. 


The  complete  set  of  substructure  gridpoint  displacement  degxees  of  freedom  is 
partitioned  to  reflect  the  division  into  essential  {*2.}  and  superfluous  |8  subsets. 
All  degrees  of  freedom  that  reside  on  interfaces  with  adjacent  substructures  must  be 
regarded  as  essential  to  the  proper  interconnection  of  substructures.  Partitioning  of  the 
displacement  set  implies  a  corresponding  partitioning  of  the  total  potential  energy  from 
Phase  2  as 
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By  definition,  the  ^  ^2bl  superfluous  to  the  objecti^re  structural  dynamics 
analyses.  This  being  the  case,  they  are  condensed  from  the  model  via  the  static  principle 
of  potential  energy^^\  This  principle  yields  a  matrix  equation  governing  static  behavior, 
..e. 


Solution  of  this  relation  for  the  superflous  degrees  of  freedom  in  terms  of  the  essential 
degrees  of  freedom  produces  a  condensing  transformation  relation  of  the  form 


where 


(28) 


(29) 


(SO) 

(31) 


Introducing  the  condensation  transformation  of  Equation  28  into  the  energy  functions  fur¬ 
nished  from  the  Phase  2  analysis  references  these  functions  to  essential  degrees  of 
freedom.  For  example,  application  to  the  strain  energy  of  Equation  22  yields 


where 


152 


The  other  energy  functions  are  similarly  transformed  to  complete  the  Phase  3  analysis. 
Generally,  the  order  of  the  matrices  can  be  substantially  reduced  by  the  introduction  of 
this  condensation  transformation. 


PHASE  4  -  MODE  SYNTHESIS 

At  the  outset  of  this  Phase  4  analysis  the  displacement  at  any  point  in  the  substructure 
is  known  in  terms  of  the  gridpoint  displacement  degrees  of  freedom  {*3}  .  These  grid- 
point  displacement  degrees-of- freedom  arise  naturally  out  of  the  finite  element  idealization 
technique  rather  than  by  choice  as  those  most  appropriate  to  structural  dynamics  analyses. 
Intuitively,  the  best  substructure  degrees -of- freedom  for  the  objective  analyses  are  those 
associated  with  the  natural  vibration  mode  shapes  of  the  substructures.  Transformation 
to  such  degrees-of- freedom  is  adopted  as  the  immediate  objective  of  this  phase. 


The  transformation  to  substructure  mode  shape  degrees-of-freedom  is  initiated  by 
partitioning  the  gridpoint  displacement  degrees-of-freedcm  into  a  subset  associated  with 
the  interface  gridpoints  ^  8  and  a  subset  associated  with  the  interior  grldpolnts 

where  the  former  are  necessarily  retained  to  effect  the  Interconnection  of  adjacent  sub¬ 
structures.  This  division  into  subsets  yields  a  corresponding  partitioning  of  the  total 

I 

potential  energy,  i.e., 


♦u4[L*3aJ’U3biJ  KJ 

Kb]' 

[[*-]•  I MJ  /M- 


Kb] 

[■'Sbb] 


(34) 


The  subset  of  displacements  at  the  interior  gridpoints  {«3b}  is  taken  to  be  made 
up  of  a  contribution  due  to  interface  displacements  |  ^  and  a  contribution  due  to  dis¬ 
placement  relative  to  the  interface 


(35) 
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A  dependence  of  interior  displacement  upon  interface  displacement  is  readily 
derived  from  static  considerations.  The  expression  of  this  dependence  follows  by  taking 
the  variation  of  the  potential  energy  (Equation  34)  and  solving  the  resulting  relation  to 
obtain 

where 

and 

{^b}  '{[sbbl'^Kb}}- 

Definition  of  displacement  relative  to  the  interface,  completes  the  expression  of 
contributions  to  the  total  displacement  at  a  point  in  the  interior  of  a  substructure.  These 
displacement  contributions  {*r}  are  constructed  of  substructure  vibration  mode  shapes. 
The  relation  governing  vibration  of  the  substructure  is  drawn  from  the  partitioned  strain 
energy  of  Equation  34  and  the  associated  kinetic  energy  to  obtain 


(36) 


(37) 


(38) 


[^3bb  ]  { ^3b}  ‘ "  [  ^3bb]  {  ^3b} ' 


(39) 


The  suppression  of  the  interface  displacements  4  S  ^  was  tacitly  assumed  in 

.  -  I  J 

the  statement  of  Equation  39.  Therefore,  the  |  8  may  be  interpreted  as  the  .  8 ^  j, 
of  Equation  35  in  this  relation.  Explicit  expression  is  given  to  ^8  ^  J  in  terms  of 
participation  factors  -f  8  \  on  vibration  mode  shapes  by  extracting  eigenvectors  from 
Equation  39  and  writing 

This  result  completes  the  development  required  to  construct  a  transformation  to 
the  final  substructure  degreea-of- freedom.  The  result  is. 


or  symbolically, 


f{°}\  ■ 

Hr 

(41) 


(42) 


Introduction  of  this  transformation  Into  the  strain  energy  form  of  Equation  21  yields 


♦u4Kj[''4]K}- 

where 


(43) 


(44) 


The  other  energy  forms  are  correspondingly  transformed.  The  results  of  this  trans¬ 
formation  yield  the  component  parts  of  a  mathematical  model  for  the  substructure  that 
is  ideally  suited  to  combination  with  other  substructures  for  structural  dynamics  analyses 
of  large  scale  systems. 

PHASE  5  -  SYSTEM  ASSEMBLY 


Mathematical  models  are  prepared  for  all  of  the  substructures  by  utilizing  the 
procedures  described  in  Phases  1  through  4.  Assembly  of  these  substructure  models 
to  form  a  model  for  the  complete  structure  is  effected  in  this  fifth  phase. 

The  assembly  of  substructures  into  a  complete  structure  proceeds  in  the  same 
manner  as  the  assembly  of  finite  elements  into  a  substructure  that  was  described  in  the 
Phase  2  analysis.  The  degrees-of-freedom  associated  with  gridpoints  common  to  more 
than  one  substructure  are  assumed  to  be  referenced  to  compatible  coordinate  axis  direc¬ 
tions.  This  being  the  case,  the  displacement  set  for  the  complete  structure  is  related  to 
that  of  its  substructure  through  a  Boolean  transformation  of  the  form 

{O' [r  5“’]  {85}-  <“=> 
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Introduction  of  this  displacement  coordinate  transformation  into  the  strain  energy  function 
of  the  substructure  as  stated  in  Equation  43  yields  this  energy  with  reference  to  the 


Summation  of  these  individual  substructure  strain  energies  yields  a  strain  energy  function 
for  the  complete  system  as 


where 


[■^5]= 


The  other  system  level  energy  functions  are  similarly  constructed  by  application  of  the 
transformations  of  Equation  45. 

PHASE  6  -  CONDENSATION  (SYSTEM  LEVEL) 

The  energj’  functions  made  available  from  Phase  5  are  expressed  in  terms  of  the 
selected  substructure  vibration  mode  shape  degrees-of- freedom  and  degrees-of- freedom 
common  to  substructure  boundaries  or  interface.  These  latter  degrees-of-freedom  were 
retained  as  essential  through  the  preceeding  analysis  phases  in  order  to  permit  the  proper 
interconnection  of  the  substructures.  The  interconnection  has  been  accomplished  at  this 
point  and  the  appropriateness  of  these  degrees-of-freedom  bears  examination. 

As  in  the  case  of  the  degrees-of-freedom  associated  with  gridpoints  interior  to  a 
substructure,  many  interface  degrees-of-freedom  will  exist  in  consequence  of  the  natural 
breakdown  of  the  structure  into  finite  elements.  In  fact,  the  degrees-of-freedom  associated 
with  these  interface  gridpoints  can  comprise  the  major  portion  of  the  complete  set.  This 
being  the  case,  it  is  worthwhile  to  condense  out  those  regarded  as  superfluous. 

Proceeding  as  in  the  substructure  related  Phase  3  condensation,  the  complete  set  of 
degrees-of-freedom  is  partitioned  into  essential  ^  superfluous  subsets. 


Projecting  this  partitioning  forward  into  the  strain  energy  of  Equation  47  and  including 
the  corresponding  external  w'ork  function  yields  the  total  potential  energy  as 

Kb] 

I 

T 


4> 


[^5bb] 


{U 

- 


(49) 


-IMNipl- 

Given  that  the  {*ab}  are  superfluous  to  structural  dynamics  analyses,  the 
stationary  principle  of  potential  energy  is  invoked  as  a  rational  means  of  establishing  a 
functional  dependence  of  these  degrees-of- freedom  upon  the  essential  degrees-of- freedom. 
The  stationary  conditions  of  the  potential  energy  follow  from  Equation  49,  i.e., 


(50) 


Solution  of  this  relation  for  the  superfluous  degrees-of-freedom  permits  construc¬ 
tion  of  the  desired  condensation  transformation  in  the  form 


I'^saa]  'M 

_Kb]"lKb]_ 

{h.} 

k=4 

— \ 

where 


(51) 

{h}-{ 

h.} 

(52) 
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[0 
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(53) 
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(54) 
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Introducing  the  transform  relation  of  Equation  51  into  the  energy  functions  furnished 
from  Phase  5  references  these  functions  to  the  essential  degrees-of-freedom  {M-  For 
example,  the  strain  energy  of  Equation  47  is  transformed  to 

♦u=iKJK]{*s} 

where 

W 

The  analysis  phases  up  to  this  point  have  taken  a  set  of  substructures  based  on 
finite  element  models  and  has  evolved  these  models  into  a  system  model  expressed  in 
terms  of  relatively  few  highly  specialized  degrees-of-freedom  considered  essential  for 
vibration  and  dynamic  response  analyses. 

PHASE  7  -  APPLICATION  OF  BOUNDARY  CONDITIONS 

Singularities  may  exist  in  the  stiffness  matrix  associated  with  the  set  of  energy 
functions  derived  in  Phase  6.  These  singularities  may  stem  from  having  deferred  the 
application  of  certain  boundary  conditions  in  order  to  broaden  the  generality  of  the  model. 
Note  that  in  Phase  2  provision  was  made  to  apply  physical  boundary  conditions  by  striking 
out  rows  and  columns  associated  with  displacements  prescribed  as  zero.  If,  however, 
the  application  of  boundary  conditions  was  vieferred,  their  introduction  can  be  accomplished 
by  a  transformation  relation  that  simply  suppresses  the  associated  gridpoint  displacement 
degrees-of-freedom.  Symbolically, 


Introduction  of  this  transformation  into  the  energy  functions  furnished  from  Phase  6, 
yields  the  following 


2  H  [“7]  {*7} 


(58) 
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where 


[s]  =  [r,]lK,][r,] 

{*,}  is  the  vector  '{8g|  with  the  bounded  degrees-of-freedom  excluded. 


PHASE  8  -  RIGID  BODY  MODE  SWEEP 


(59) 


The  singularities  that  are  more  difficult  to  deal  with  are  those  associated  with 
unrestrained  response  of  the  structure.  These  are  accounted  for  in  the  following  para¬ 
graphs.  First,  a  set  of  rigid  body  modes  is  constructed.  Subsequent  simplification  is 
realized  if  these  are  referenced  to  the  center  of  gravity  of  the  structure  but  this  is  not 
essential.  The  statement  of  a  set  of  rigid  body  modes  is  particularly  simple  for  the  type 
of  model  developed  because  none  of  the  degrees-of-freedom  which  are  amplitude  coefficients 
of  substructure  vibration  mode  shapes  participate  in  a  rigid  body  motion.  Only  the  rela¬ 
tively  few  gridpoint  displacement  degrees-of-freedom  that  have  been  retained  at  the  sub¬ 
structure  interfaces  need  be  considered.  The  desired  rigid  body  modes  are  easily  formed 
and  are  indicated  symbolically  herein  as  [  R  3  . 

Excluding  damping  and  excitation  forces  leaves  only  strain  and  kinetic  energies 
from  which  the  matrix  equation  of  motion  can  be  readily  extracted,  i.e., 

[m,]{8,}  +  [k,]{8,}  =  {o}.  (60) 

Premultiplication  of  this  result  by  the  transpose  of  the  rigid  body  mode  matrix  [R] 

yields 


This  multiplication  negates  the  second  term  and  permits  direct  integration  to  obtain  a 
relation  governing  rigid  body  motion  as 


(62) 
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The  rigid  body  motion  of  Equation  62  may  be  neglected  (|oi}  =  { }  “  {oj)  for 
present  purposes  and  the  coefficient  matrix  written  In  partitioned  form  using  a  new 
symbol  for  convenience,  l.e., 


(63) 


The  selection  of  the  degrees-of-freedom  for  Inclusion  Into  |8  Is  largely  a 
matter  of  convenience  within  the  constraint  that  the  [  be  square  and  nonsingular. 

The  order  of  y  j  Is  equal  to  the  number  of  rigid  body  modes.  Solution  of  Equation  63 
for  the  •  In  terms  of  the  permits  construction  of  the  transformation  sought 

to  extract  the  rigid  body  modes  from  the  model  furnished  from  Phase  7,  l.e.. 


where 


{^s}"  t^7a} 


1  - 

[o] 

• 

J 

[-[ 

1" 

[^a] 

(64) 


(65) 


This  Phase  8' of  the  analysis  process  Is  completed  by  the  Introduction  of  the  trans¬ 
formation  [  defined  In  Equation  65, 

The  result  for  the  strain  energy  function  Is 


where 

K]  =  [^rK]W- 


(66) 


(67) 


All  other  energy  functions  transform  correspondingly.  Note  that  this  phase  can  be 
neglected  if  the  physical  boundary  conditions  in  Phase  7  are  such  as  to  prevent  rigid 
body  motion. 
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PHASE  9  -  SYSTEM  VIBRATION  MODES  AND  FREQUENCIES 

The  mathematical  model  has  been  brought  to  the  point  where  undaniped  vibration 
modes  and  frequencies  can  be  determined.  The  order  of  this  eigenvalue  problem  is 
presumed  to  have  been  reduced  to  well  within  program  capacity.  The  relevant  governing 
relation,  derivable  from  the  strain  energy  and  kinetic  energy,  is  given  by 

This  relation  is  rewritten  in  the  form 


(69) 


where 


(70) 


to  facilitate  extraction  of  the  mode  shapes  corresponding  to  the  lowest  natural  frequencies. 
These  natural  frequencies  and  mode  shapes  of  vibration  complete  the  characterization 
required  by  some  design  specifications. 


This  phase  of  the  analysis  is  extended  here  to  carry  the  computations  forward  to 
provide  a  basis  for  prediction  of  time  dependent  response.  The  eigenvectors  of  Phase  8 
are  collected  together  to  form  the  columns  of  a  modal  matrix  designated  [  Tg  j .  Follow¬ 
ing  the  normal  mode  approach  to  dynamics  analysis,  ( )  this  matrix  is  employed  to  effect 
a  transformation  to  degrees-of-freed  jm  M  which  are  participation  coefficients  of  the 
natural  mode  shapes  of  the  complete  structure,  i.e., 

{«8}=[r8]K}-  <”) 


This  transformation  completes  the  development  of  an  optimum  set  of  degrees-of- 
freedom  for  use  in  the  prediction  of  dynamic  response.  The  associated  final  form  of  the 
strain  energy  is  derived  by  substitution  into  Equation  66  to  obtain 

[^9]  [^9]  {^9} 
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where 


K]=  [r9]''[''8][rj.  P3) 

At  this  point  it  is  instructive  to  reconstruct  the  complete  sequence  of  transforma¬ 
tions  that  are  applied  to  a  typical  finite  element  function  to  arrive  at  this  final  form. 

Working  backwards,  the  final  system  level  quadratic  forms  stem  from  a  three  stage 
transformation  beyond  the  point  of  assembly  of  the  substructures.  In  order,  these  are 
condensation  [  Tg  j  (Equation  53),  singularity  sweep  [  Tg  j  (Equation  65)  and  normal 
mode  [  ^g]  (Equation  71)  transformations.  These  yield  a  collective  modification  of  the 
system  level  stiffness  matrix  in  going  from  the  /  to  the  that  is  given  by 

[K9]  =  [r9]Mr3]nr,]nre]^k][re][r,][r9][r9] 

The  contribution  of  each  substructure  to  the  initial  system  level  stiffness  matrix 
[  ^5^^^  ]  derived  from  transformations  beyond  the  point  of  assembly  of  its 

component  finite  elements.  In  order,  these  are  prescribed  displacement  ^  F  ^  j  (Equation 
20),  condensation  l^r  gj  (Equation  30),  component  modes  [F^j  (Equation  42)  and 
assembly  [F  g  j  (Equation  45)  transformations.  These  transformations  yield  a  collective 
modification  of  the  substructure  level  stiffness  matrix  in  going  from  the  to  the 

that  is  given  by 

[■'a]  =  I  [r3]lr,]-[K,][r,]  [r,]  [r,]  [r,]. .. 

Tracing  the  sequence  of  transformations  to  the  fundamental  finite  element  blocks  is 
completed  by  statement  of  the  element  assembly  transformation  [  F  j  of  Equation  13,  i.e. 

[k®]  =  Z  [r«]^[K«][r<'’>]  (76, 

This  recap  of  the  dynamic  substructuring  procedure  makes  clear  the  considerable 
computation  involved.  It  also  exhibits  the  highly  systematic  nature  of  the  process. 


;4<^jW!!m^lJ^4^^J^^^■J;^V.^!'^WiAU»j^'■^u  I.JMHUI.,'  .  "t  "'- 


>W»f»^  '**}■'**»  H<I»WK»» 


NOTATION 

Rectangular  matrix  (£q  1) 

Column  matrix  (Eq  1) 

Row  matrix  (Eq  2) 

Summation  operator  (Eq  16) 

Finite  element  displacement  functions  (Eq  1) 

Finite  element  assumed  displacement  mode  shapes  (Eq  1) 

Gridpoint  displacement  coefficients  of  finite  element  assumed 
displacenient  mode  shapes  (Eq  1) 

Strain  energy  function  (Eq  2) 

External  work  function  (Eq  3) 

Total  potential  energy  function 

Finite  element  stiffness  matrix  (Eq  2) 

Finite  element  applied  load  vector  (Eq  3) 

Psuedo  potential  for  viscous  damping  (Eq  4) 

Psuedo  potential  for  structural  damping  (Eq  5) 

Unit  imaginary  number  (Eq  7) 

Finite  element  viscous  damping  matrix  (Eq  4) 

Finite  element  structural  damping  matrix  (Eq  5) 

Finite  element  kinetic  energy  function  (Eq  6) 

Differentiation  with  respect  to  time  (Eq  6) 

Mass  matrix  (Eq  6) 

Gridpoint  displacement  coefficients  referenced  to  a  coordinate 
system  defined  by  the  subscript  £  (Eq  8) 

Gridpoint  displacement  coefficients  referenced  to  a  coordinate 
system  defined  by  the  subscript  ^+1  (Eq  8) 

Translational  transformation  which  relates  gridpoint  displacement 
coefficients  referenced  to  a  coordinate  system  defined  by  the 
subscript  £  to  one  defined  by  -/+  1 
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[r/.  l] 


Rotational  transformation  which  relates  gridpoint  displacement 
coefficients  referenced  to  a  coordinate  system  defined  by  the 
subscript  M.  to  one  defined  by  ^  +  1 

Substructure  identification  number  (Eq  13) 


k  Finite  element  identification  number  within  a  substructure  (Eq  1,13) 


T  Matrix  transpose  operator  (Eq  10) 


th 

The  complete  set  of  gridpoint  displacements  of  the  j  substructure 
(fiq  13) 


Transformation  relating  r  to  the  displacement  set  of  the  k 

finite  element  (Eq  13)  ^ 

The  subset  of  substructure  displacements  that  are  prescribed  (Eq  17) 

The  subset  of  substructure  displacements  that  are  essential  degrees 
of  freedom  (Eq  17) 


[''3] 

{*3.} 

{^3.} 


A  contribution  to  the  external  work  which  arises  from  the  prescribed 
displacements  (Eq  24) 

The  subset  of  substructure  gridpoint  displacement  degrees  of  freedom 
that  is  regarded  as  essential  to  structural  dynamic  aralyses  (Eq  26) 


The  subset  of  substructure  gridpoint  displacement  degrees  of  freedom 
that  is  regarded  as  superfluous  to  structural  dynamics  analysis  (Eq  26) 


Generalized  loads  corresponding  to  the  <  $  > 

Generalized  loads  corresponding  to  the  ,  . 


(Eq  26) 
(Eq  26) 


Rotational  transfcrmatlx^n  relating  the  complete  set  of  substructure 
gridpoint  displacement  degrees  of  freedom  to  the  subset  essential 
to  structural  dynamics  analyses  (Eq  28) 


Translational  transformation  relating  the  complete  set  of  substructure 
gridpoint  displacement  degrees  of  freedom  to  the  subset  essential  to 
dynamics  analyses  (Eq  28) 


Substructure  stiffness  matrix  referenced  to  the  gridpoint  displacement 
degrees  of  freedom  essential  to  the  structural  dynamics  analyses 
(Eq  33) 

{^3}  displacements  associated  with  gridpoints  interior  to  a 

subsiructure  (Eq  34) 

5^3^  subset  of  displacements  associated  with  gridpoints  on  the  interface 
a  pv*bstructure  (Eq  34) 
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Contribution  to  {«  3b  }  •rising  from  displacements  of  the  interface 
gridpoints  of  a  substructure  (£q  351 

Contribution  to  {»  3b  }  arising  from  displacements  relative  to  interface 
gridpoints  of  a  substructure  (Eq  35) 

Transformation  relating  interface  displacements  to  the  displacements 
which  these  induce  at  (Eq  37)  interior  points. 

Displacement  contribution  at  interior  points  due- to  implied  loads  at 
interior  points  (Eq  38) 

Modal  matrix  of  substructure  taken  as  constrained  at  the  substructure 
interfaces  (Eq  40) 

Participation  jfactors  on  substructure  mode  shapes 

Rotational  tiansformation  from  essential  gridpoint  degrees  of  freedom 
to  the  final  set  of  substructure  degrees  of  freedom  (Eq  42) 

Translational  transformation  from  essential  gridpoint  degrees  of  freecc  m 
to  the  final  set  of  substructure,  degrees  of  freedom  (Eq  42) 

Final  set  of  substructure  degrees  of  freedom  (Eq  42) 

Final  form  of  substructure  stiffness  matrix  for  assembly  into  complete 
structure  (Eq  44) 

Complete  set  of  degrees  of  freedom  after  assembly  of  the  substructures 
(Eq  45) 

Transformation  between  final  substructure  degrees  of  freedom  and  the 
initial  set  of  degrees  of  freedom  for  the  assembled  structure  (Eq  45)  ' 

Stiffness  matrix  of  complete  structure  referenced  to  the  Js  I 
(Eq  48) 

The  subset  of  degrees  of  freedom  in  |Scl  regarded  as 

essential  (Eq  49)  • 

The  subset  of  degrees  of  freedom  in  is  regarded  as 

siqperfluous  (Eq  49)  J 

The  applied  load  vector  corresponding  to  the  degrees  of  freedom  <  8.  ^ 
(Eq  49).  I  5aJ 

The  )9^olied  load  vector  corresponding  to  the  degrees  of  freedom  /  8_,  I 
,(Eq,49) 

Condensation  transformation  from  the  ^8  gj  degrees  of  freedom  for  the 
total  structure  to  the  subset  chosen  as  essential  (Eq  53) , 

Contributions  to  superfluous  degrees  of  freedom  from  the  corresponding 
applied  loads  (Eq  54)  ' 
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SECTION  V 


DISCUSSION  AND  CONCLUSIONS 


A.  DISCUSSION 

Integrated  general  purpose  analysis  capabilities  of  the  MAGIC  II 
System  (ilass  signal  a  major  advance  in  the  state-of-the-art  of  automated 
tools  for  analysis.  The  superior  cost  effectiveness  of  such  systems 
over  conventional  multiple  special  purpose  program  capabilities  is 
compelling. 

This  assertion  of  superior  performance  from  large  scale  program 
systems  may  well  contradict  conclusions  drawn  from  experience. 

Complexity  and  inefficiency  have  long  been  concomitant  with  large 
size  and  versatility  iii  computer  programs.  Indeed,  the  elimination 
of  these  depreciating  effects  was  prerequisite  to  realization  of 
the  favorable  cost  effectiveness  of  the  MAGIC  II  System. 

Large  size  and  versatility,  without  excessive  complexity,  are 
assumed  Intrinsic  to  the  MAGIC  II  System  in  subsequent  paragraphs,  as 
attention  is  focused  upon  the  relative  efficiencies  of  integrated 
general  purpose  analysis  capabilities  and  multiple  special  purpose 
computer  program  analysis  capabilities.  This  is  to  presume  the 
pre-requisite  elimination  of  the  greater  hindrance;  namely,  the 
excessive  complexity  which  choked  off  many  early  general  purpose 
program  developments.  This  problematical  complexity  was  encountered 
when  programs  of  simple  organization  grew  to  press  upon  the  limits 
of  computer  software  £ind  hardware  capabilities.  Extensions  beyong 
this  point  were  accomplished  by  intricately  coordinated  multiple  usage 
of  valuable  names  and  locations,  special  program  versions  with  omitted 
features  and  other  actions  which  accumulated  to  entangle  the  logic  and 
data  storage  until  further  modification  became  Impractical. 

In  the  face  of  this  situation  increasingly  powerful  analytical 
models  and  solution  methods  were  formulated  and  numerical  implemen¬ 
tation  demanded.  And,  as  is  often  the  case,  sufficient  pressure  was 
built  up  to  bring  about  the  technological  advances  needed  in  the 
computer  technologies . 


i^dvances  were  forthcoming  in  programming  technology  which 
established  the  technical  feasibility  of  a  truly  general  purpose 
computer  program  system.  Advances  in  computer  hardware  insured  the 
economic  feasibility  as  the  technical  feasibility  was  established 
through  a  number  of  contributing  developments.  The  collective 
result  of  these  latter  developments  is,  in  a  word,  "organization". 

Among  those  organizational  characteristics  or  features  considered 
essential  are,  the  breakdown  into  single  function  modules,  the 
program  library  concept,  the  matrix  interpretive  system,  machine 
independency,  etc. 

It  is  appropriate  to  emphasize  at  this  point,  that  the 
MAGIC  II  System  for  structural  analysis  is  more  than  a  discrete 
element  computer  program.  It  is,  in  one  sense,  a  Problem  Oriented 
Language  (POL)  which  enables  various  Analyst  specified  computa¬ 
tional  procedures.  And,  at  the  same  time,  it  is  designed  with  attendant 
structural  analysis  practices  evolved  from  applications  experience. 

These  practices  are  discussed  in  detail  in  subsequent  paragraphs. 

The  point  of  interest  here  is  that  the  efficiency  of  the  MAGIC  II 
System  is  an  overall  efficiency  governed  more  by  men  than  machines. 

The  more  comprehensive  the  comparison,  the  greater  the 
advantage  shown  by  the  integrated  general  purpose  analysis  capabili¬ 
ties  over  multiple  special  purpose  program  capabilitites.  In  nearly 
all  cases  an  equitable  comparison  must  include  consideration  of 
program  development  efforts  since  relevant  technologies  are 
continuously  advanced.  On  this  basis  the  integrated  approach 
enjoys  the  greatest  relative  advantage.  The  integrated  approach  is 
also  superior  to  the  multiple  program  approach  when  considering  only 
factors  involved  in  utilization  of  operational  capability.  On  the 
other  hand,  shorter  execution  times  are  conceded  to  special  purpose 
programs  without  dispute,  since  execution  efficiency  is  not  essential 
to  the  case  for  the  greater  overall  efficiency  of  integrated  analysis 
capabilities . 


Attention  is  focused  now  on  the  impact  of  the  integrated 
general  purpose  computer  program  approach  on  the  efficiency  of 
the  many  processes  involved  in  maintenance  and  application  of 
responsive  analysis  tools  in  support  of  a  broad  structural  design 
activity.  Program  maintenance  efforts  benefit  from  the  highly 
modularized  organizational  structure  to  an  even  greater  extent 
than  the  initial  development  effort. 

In  the  Initial  development,  functional  modules  are 
established  against  the  requirements  of  the  alternative  analysis 
procedures  taken  collectively.  And,  since  an  extensive  commonality 
exists,  multiple  repetitious  coding  is  avoided.  This  same 
payoff  is  derived  again  as  existing  modules  are  retired  in  favor 
of  new  modules  which  offer  improved  performance.  The  introduction 
of  a  single  improved  module  is  reflected  to  advantage  throughout 
all  pertinent  analysis  procedures  of  the  computer  program  system. 

The  option  exists  to  retain  alternative  modules  for  the  same 
function  without  sacrifice.  This  provides  useful  operational 
flexibility  and  a  convenient  testbed  for  various  candidate  procedures. 
Alternative  procedures  can  be  evaluated  within  the  system  without 
disrupting  its  operational  status. 

The  foregoing  has  dealt  with  maintenance  of  existing  analysis 
capability.  Maintenance  is  also  interpretable  as  generalization  of, 
and  addition  to,  the  overall  analysis  capability.  Completely  new 
analyses  can  be  implemented  with  the  addition  of  only  those  functional 
modules  absent  in  the  existing  capability.  For  example,  finite 
element  heat  conduction  analyses  and  the  more  efficient  optimality 
criteria  based  optimization  methods  are  possible  with  relatively 
minor  modifications  to  the  f4AGIC  II  System. 

The  benefits  derived  from  the  organization  of  a  general  purpose 
computer  program  system  in  development, .maintenance,  generalization 
and  extension  are  simultaneously  important  disadvantages  associated 
with  multiple  computer  program  analysis  capabilities.  The  extensive 
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commonality  among  analyses  leads  in  this  latter  case  to  the 
repeated  development  of  coding  to  perform  a  given  function.  The 
preparation  of  special  versions  of  new  modules  and  the  introduction  of 
these  into  a  multiplicity  of  computer  programs  is  often  not  justified 
and  the  overall  capability  is  depreciated. 

Another  particularly  Important  handicap  borne  by  the  separate 
programs  of  a  multiprogram  capability  is  that  these  programs  cannot 
command,  individually,  the  provision  of  many  useful  special  features. 
For  example,  useful  options  and  diagnostics  are  usually  omitted  from 
these  special  purpose  program  routines.  Also,  such  programs  frequently 
encounter  obstacles  such  as  machine  storage  capacity  which  must  be 
avoided  rather  than  surmounted  in  view  of  the  limited  applicability 
of  the  program.  Advancements  in  computer  Sf^ftware  and  hardware  are 
further  considerations  of  importance  in  the  maintenance  of  an  analysis 
capability.  These  advancements  place  multiple  program  capabilities 
in  special  peril.  Those  programs  not  being  actively  utilized  at  the 
time  of  transition  in  software  or  hardware  are  easily  overlooked  and 
in  this  way  are  lost  from  the  overall  analysis  capability. 

No  single  factor  is  more  important  in  the  provision  of  a 
responsive  analysis  capability  than  documentation.  Engineering 
documentation  must  delineate  analysis  procedure,  input  data  and 
output  data.  Programming  documentation  must  provide  for  operation 
and  modification  of  the  program. 

Consolidation  of  the  analysis  capability  into  a  general  purpose 
program  results  in  a  corresponding  favorable  consolidation  of 
documentation.  Not  only  is  volume  reduced  but  the  total  capability 
is  described  uniformly  as  a  whole.  Small  programs  tend  to  be  the 
personal  tool  of  the  initiator.  As  a  consequence,  the  documentation 
prepared  is  generally  inadequate  to  enable  general  usage.  This 
situation  leads  to  extensive  tutorial  instruction  to  realize  the 
benefits  of  the  program  development.  At  the  very  least,  multiple 
program  capabilities  place  the  burden  of  assimilating  rhe  overall 
analysis  capability  from  the  individual  manuals  upon  tl-e  user. 
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The  foregoing  has  pointed  out  decisive  advantages  of  general 
purpose  program  systems  in  the  context  of  development  and 
maintenance  of  analysis  capability.  The  most  compelling 
advantages,  however,  are  found  in  operation.  The  greater  efficiency 
of  the  MAGIC  II  System  relative  to  multiprogram  capabilities  for 
analysis  stems  in  large  measure  from  the  extent  of  the  analysis 
process  which  is  covered.  Time  consuming,  error  prone,  manual 
transfers  of  data  between  special  purpose  or  single  step  computer 
programs  are  avoided.  The  Integration  of  heat  conduction  and 
thermal  stress  analysis  within  a  single  system  can  circumvent  the 
laborious  preparation  of  temperature  data.  The  integration  of 
stiffness  and  vibration  analyses  can  similarly  circumvent  the 
manual  transfer  of  stiffness  and  mass  data.  These  eliminations  of 
manual  effort  yield  reductions  in  calendar  time  which  is  often  the 
paramount  consideration  for  contribution  of  analysis  to  design. 

This  is  not  to  say  that  long  continuous  executions  are  desirable. 
Execution  interruptions  enter  importantly  into  proper  utilization 
of  the  MAGIC  II  System. 

The  MAGIC  II  System  is  designed  to  facilitate  good  structural 
analysis  practices  in  support  of  the  overall  structural  design 
process.  Individual  design  orgaiiizations  are  best  served  by 
structural  analysis  practices  and  program  versions  which  are,  to 
some  degree,  distinct.  On  the  other  hand,  the  extensive  commonality 
which  does  exist  among  design  organizations  provides  strong 
motivation  for  reviewing  the  effective  structural  analysis  practices 
and  supplemented  program  version  which  have  evolved  at  Bell  Aerospace 
Company. 

The  structural  analysis  process  begins  with  the  idealization  of 
the  structure  into  an  assemblage  of  finite  elements.  This  is  a 
multistep  operation  if  the  structure  is  first  separated  into  sub¬ 
structures.  Generally,  the  separation  into  substructures  is 
governed  by  the  physical  interconnections  of  the  major  structural 
components.  The  idealization  into  finite  elements  is  governed  by 
variations  in  geometry,  dimensions,  material,  applied  loading  and 
boundary  conditions. 
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Preprinted  input  data  forms  are  employed  to  simplify  and  thereby 
improve  the  reliability  of  the  input  data  specification.  These 
preprinted  input  forms  associated  with  the  MAGIC  II  System  are  an 
important  improvement  over  card  image  forms  for  frequent  as  well  as 
Infrequent  users  since  they  incorporate  automatic  data  generation 
features.  These  built-in  data  generation  features  are  supplemented 
at  Bell  by  auxiliary  (not  Integrated  into  the  MAGIC  II  System)  data 
generation  programs.  Some  of  these  are  employed  routinely.  Others 
are  extremely  simple  programs  written  for  a  single,  problem 
related  calculation.  Such  auxiliary  programs  are  frequently  employed 
to  advantage  in  the  generation  of  grldpoint  coordinates  with 
reference  to  the  global  rectangular  coordinate  axes,  since  expression 
of  these  can  require  extensive  tedious  calculation.  This  gridpoint 
coordinate  data  set  should  be  interpreted  here  to  include  points 
for  specification  of  grldpoint  axes  transformations  and  stress  and 
material  angles  as  well  as  points  associated  with  degrees -of - 
freedom. 

The  first  MAGIC  II  System  execution  undertaken  is  to  confirm 
the  assembled  input  data  deck.  This  deck  is  read  and  the  implied 
data  is  given  explicit  definition.  For  example,  material  properties 
are  extracted  from  the  Material  Library  and  gridpoint  axes  trans¬ 
formations  are  generated  from  the  coordinate  table.  The  completed 
data  set  is  examined  in  this  preprocessing  execution.  All  data 
items  are  stored  for  execution  restart  and  printed  for  further 
checking  by  the  analyst. 

The  preprocessing  execution  is  supplemented  at  Bell  to  include 
the  generation  of  a  magnetic  tape  which,  in  turn,  generates  a  plot 
of  the  structural  model  automatically.  This  plot  enables  efficient 
and  reliable  confirmation  of  the  two  most  problematical  data  items j 
namely,  the  gridpoint  positions  and  the  finite  element  connection 
arrangement.  Beyond  this  point  the  structure  plot  is  a  useful 
identifying  title  sheet  for  the  printed  problem  output. 
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The  next  phase  of  the  analysis  process  proceeds  via  a  restart 
through  the  generation  of  the  structural  matrices  for  stiffness, 
stress,  loads,  assembly,  boundary  conditions,  etc.  Built-in 
features  control  this  matrix  generation  to  selectively  form  only 
those  matrices  required  for  the  current  analysis.  Completion  of 
the  matrix  generation  phase  signals  exit  from  the  Structural  System 
Monitor.  This  is  an  interface  point  between  matrix  abstraction 
Instruction  statement,  and,  therefore,  a  point  for  optional 
interruption  of  the  execution  to  examine  the  system  level  matrices. 
This  interruption  is  used  only  infrequently  at  Bell. 

Calculation  proceeds  to  the  governing  matrix  equation  and 
thence  to  the  solution  for  the  displacement  vectors  for  all  load 
conditions.  For  some  problems  execution  may  be  terminated  at  this 
point.  For  many  other  problems  the  validity  of  the  analysis  can  be 
assessed  against  these  displacement  results  and  an  execution  interrup¬ 
tion  is  justified  by  the  computational  investment  required  for  the 
secondary  results,  ideally,  the  deformed  structure  should  be  plotted 
to  facilitate  interpretation  of  the  predicted  displacement  behavior. 

The  analysis  proceeds  from  the  displacement  solution,  with  or 
without  interruption,  to  calculation  and  print  of  the  remainder  of 
the  output  data  items j  namely,  reactions,  forces,  stresses,  etc.  This 
is  the  conventional  point  of  termination  of  finite  element  analyses. 
However,  a  number  of  relatively  simple  auxiliary  programs  are  used 
to  advantage  at  Bell  to  relieve  the  burden  this  output  places  on  the 
stress  engineers.  As  in  the  case  of  the  input  data  generation 
auxiliary  programs,  some  of  the  auxiliary  output  data  reduction 
programs  are  employed  repeatedly  and  others  are  special  to  a  single 
problem.  The  functions  of  these  programs  Include  such  things  as 
principal  stress  calculations  and  margin  of  safety  determinations. 
Auxiliary  programs  which  do  nothing  but  selectively  print  and  label 
output  data  items  are  also  helpful  for  large  problems. 


Several  comments  on  the  evaluation  of  output  data  are 
warranted  in  concluding  discussion  of  good  structural  analysis 
practices.  The  examination  of  output  by  the  Analyst  should  be 
initiated  under  the  presumption  that  an  error  exists  with  confidence 
in  the  validity  of  the  analysis  accumulating  as  the  examination 
proceeds.  Given  a  complete  set  of  output,  attention  should  first 
be  given  to  the  gridpoint  force  balances  and  reactions.  Assured  that 
no  unintended  reactions  exist  and  that  residuals  are  negligibly 
small,  the  displacement  states  should  be  examined.  If  the  general 
deformed  configuration  does  not  expose  any  inconsistencies, 
confirmation  is  completed  by  examination  of  the  more  extensive 
presentation  of  force  and  stress  data. 

The  foregoing  discussion  has  focused  upon  development, 
maintenance  and  utilization  considerations  important  to  the 
favorable  cost  effectiveness  of  the  present  MAGIC  II  System 
for  structural  analysis.  Further  evolution  of  this  system 
can  be  expected  which  will  continue  to  improve  its  relative 
advantage.  Updated  versions  of  the  MAGIC  II  Sys cem  will  be 
compatible  with  all  features  developed  in  connection  with  prior 
versions. 

B.  CONCLUSIONS 

It  is  concluded  that  the  fviAGIC  II  System  is  a  logical  and 
consistent  extension  of  the  original  MAGIC  System  and  that 
additional  capabilities  realized  ;*»ith  the  System  have  met  or 
exceeded  the  requirements  of  Contract  F  336l5-69-C-124l.  The 
satisfactory  achievement  of  the  overall  objectives  is  given 
substantiation  by  a  number  of  subsidiary  conclusions.  Specifically, 
it  is  concluded  that : 

(1)  The  versatile  finite  element  library  enables  effective 
idealization  of  most  linear  structures. 

(2)  Computational  procedures  attendant  to  the  MAGIC  II  System 
enable  the  conduct  of  linear  displacement  and  stress 
'Analyses  in  the  presence  of  general  prestrain  and  thermal 
loading  as  well  as  distributed  and  concentrated 
mechanical  loading.  Additionally,  vibration  analyses  can 

be  employed  with  or  without  the  use  of  condensation  techniques. 
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(3)  The  stability  analysis  procedure  provided  in  the  MAGIC  II 
System  enables  the  prediction  of  critical  load  levels  for 
general  built-up  shell  structures. 

(4)  The  preprinted  input  data  forms  facilitate  the  rapid  and 
reliable  specification  of  problem  data  as  evidenced  by 
their  wide  acceptance  with  the  original  MAGIC  System. 

(5)  The  output  provided  by  the  MAGIC  II  System  is  oriented 

to  the  engineering  user  and  facilitates  clear  and  concise 
interpretation  of  output  parameters. 

(6)  The  computer  program  organization  of  the  MAGIC  II  System 
is  logical  in  design  and  is  wel3  suited  to  generalization. 
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